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EXPERIMENTAL AND THEORETICAL INVESTIGATION OF PASSIVE 
DAMPING CONCEPTS FOR MEMBER FORCED AND FREE VIBRATIONS 

By 

Zia Razzaq 1 and David W. Mykins 2 
ABSTRACT 

The results presented in this research report are the outcome of an 
ongoing study directed toward the identification of potential passive damp- 
ing concepts for use in space structures. The effectiveness of copper 
brush, wool swab, and n silly putty 11 in chamber dampers is investigated 
through natural vibration tests on a tubular aluminum member. The member 
ends have zero translation and possess partial rotational restraints. The 
silly putty in chamber dampers provide the maximum passive damping efficien- 
cy. Forced vibration tests are then conducted with one, two, and three 
silly putty in chamber dampers. Owing to the limitation of the vibrator 
used, the performance of these dampers could not be evaluated experimentally 
until tlie forcing function was disengaged. Nevertheless, their performance 
is evaluated through a forced dynamic finite element analysis conducted as a 
part of this investigation. The theoretical results are based on experi- 
mentally obtained damping ratios indicate that the passive dampers are 
considerably more effective under member natural vibration than during 
forced vibration. Also, the maximum damping under forced vibration occurs 
at or near resonance. 

* Professor, Department of Civil Engineering, Old Dominion University, 
Norfolk, Virginia 23529. 

^Graduate Research Assistant, Department of Civil Engineering, Old Dominion 
University, Norfolk, Virginia 23529. 
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NOMENCLATURE 


[C] = damping matrix for member 

[D] = displacement vector 

* 

[D] = velocity vector 

• • 

[D] = acceleration vector 

[Dj] =» displacement vector at node j 

E = Young’s Modulus 

I = moment of inertia 

[K] « global stiffness matrix for member 
[R] — forcing function vector 

y 6 “ arbitrary constants for Newmark’s method 

A* d = dynamic deflection amplitude 

A p = static midspan deflection 

A§ = static midspan deflection 

At « time increment 

$ =* modal vector 

n = damping efficiency index 

Q = frequency of applied forcing function 

co - natural frequency 

co fe = natural frequency from finite element analysis 
p = mass density 
£ « damping ratio 
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1 . INTRODUCTION 


1.1 Background and Previous Work 

The space station designs currently under consideration by NASA are 
three-dimensional space structures composed of long tubular members. 

Modules providing the required living and working space for astronauts will 
be attached to this framework. Such a structure, suspended in a weightless 
environment, would be subjected to many types of dynamic loading. These 
include differential heating or cooling of the structure, variations in 
acceleration or gravitational pull, and impact with a solid object. The 
ability to expeditiously damp these vibrations before they cause permanent 
damage is a practical problem worth studying. 

The necessarily large slenderness ratio of the average space truss 
member, combined with the flexible, semi-rigid end restraints cause the 
dynamic response of these members to be characterized by low frequency, 
small amplitude vibrations. Active damping techniques utilize electronic 
sensors and movable masses to reduce vibration of structures. This system, 
although effective, requires regular maintenance and an external power 
source. An alternative for mechanically damping a system is the concept of 
"passive" damping. This method uses a device or material permanently 
attached to the structure or its components and designed to absorb the 
energy of vibration thereby providing some damping of the system. Unlike 
"active" damping, this would require minimal maintenance and no external 
power. 

The challenge to developing a passive damping concept, particularly for 
a space structure is two- fold. First is the necessity to minimize the mass, 
for without this constraint one obvious solution would be to provide large 
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mass concentrations at the critical nodal points for the vibration modes. 
Such an approach would be expensive since the cost of transporting the 
system into space is directly related to the mass. The second challenge is 
to identify a concept which will provide passive damping without altering 
the strength or stiffness of the structure. For example, mild compression 
of the members provides some damping, however, the safe service loads for 
the structure are altered. 

Recently, investigations into passive damping concepts for slender 
tubular members have been conducted with various end conditions (References 
1-5) . The most effective concepts found were the mass-string-whiskers 
assembly, and brushes for electrostatic and frictional damping. In these 
experiments, only natural flexural vibration was examined. 

The previous work was conducted on hollow tubular steel members with an 
outer diameter of 0.5 inches. The passive damping concepts which were found 
to be effective for these members may not be as effective if the dimensions 
are changed. Factors altered by dimens ionsal changes may include the damper 
mass required, the extent of the frictional interaction, and the member 
dynamic characteristics. Clearly research is needed to identify a viable 
passive damping concepts for members of different sizes and dymanic 
properties. In the present study, hollow tubular aluminum members with an 
outside diameter of 2.0 inches are used. These members more closely 
resemble the actual size and material which may be used in the future space 
stations. 

1.2 Problem Definition 

Figure 1 shows schematically a slender beam of length L with a hollow 
circular cross section. The outer diameter is D 0 the inner diameter is , 
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and the material is aluminum with a Young’s modulus of 10,000 ksi. An 
aluminum member is used because the graphite composite tubes which may 
possibly be used in space structures are not yet available. The member ends 
are provided with a prototype connection developed by NASA for the space 
frames. These connections possess partial rotational restraint 
characteristics in the plane of motion and a more rigid end condition in the 
orthoganal plane. No axial or lateral movement of the member ends is 
permitted. 

The problem is to identify potential passive damping concepts to absorb 
the energy of both natural flexural vibration, and harmonic forced flexural 
vibration, and to study the effectiveness of each concept. The natural 
vibration is caused by the sudden release of a constant static load. The 
harmonic forcing function is applied through a mechanical connection to a 
harmonic vibrator. 

1.3 Objective and Scope 

The following are the main objectives of this study: 

1. Identification of potential passive damping concepts for slender 
tubular structural members. Specifically, the following damping 
concepts are investigated: 

a. wool swabs, 

b . copper brushes , 

c. silly putty in chambers. 

2. Evaluation of the damping efficiencies of the various damping 
concepts . 

3. Evaluation of the suitability of a theoretical finite element 
analysis by comparison to experimental results for natural and 
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forced vibration, and a previous finite-difference solution for 
natural vibration. 

Only flexural member vibration is considered. The natural vibration 
study is conducted on each of the three passive damping concepts and for one 
specific initial deflection. Only the most efficient damping concept is 
considered for further study under forced vibration. Also, the vibration is 
induced by load application at the member midspan. 

1.4 Assumptions and Conditions 

The following assumptions and conditions have been adopted in this 
study: 

1. The deflections are small. 

2. The material of the member is linearly elastic. 

3. Only planar vibration is considered. 

4. Damping force is opposite but proportional to the velocity at any 
location along the member. 

5. The damping force is uniform along the length of the member. 

6. The member is tested at 1-g and room temperature. 

7. The effect of secondary induced forces such as varying axial 
tension and compression developed in the member during vibration 
is considered to be negligible. 
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2 . THEORETICAL FORMULATION 


2.1 Finite Element Formulation 

The beam shown in Figure 1 may be divided into N finite elements along 
the length. For the discretized system, the governing equation of motion 
can be expressed in the following matrix form (Reference 6) : 


[K] {D} + [M]{D) + [C] (D) - {RJ (1) 

where : 

{ D ) = displacement vector, 

{ D } - velocity vector, 

{ D } — acceleration vector, 

[K] -global stiffness matrix for the "structure", 

[M] - global modified lumped mass matrix, 

[C] - damping matrix, 

{R} = forcing function vector. 


The boundary and initial conditions for the problem shown in Figure 1 
are given in Reference 2 and are summarized here: 


D(0 , t) - 0 
D(L,t) - 0 

El D" (0 , t) - D’ (0,t) 
El D" (L , t) - -k 2 D'(L.t) 
D(x,0) - 0 

D(x,0) - 0(x; K, El, L) 


( 2 ) 

(3) 

(4) 

(5) 

( 6 ) 
(7) 
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where primes represent differentiation relative to x, and dots represent 
time differentiation. The displacement vector at any node j along the 
member can be written as : 




{Dj } - 



( 8 ) 


where dj and d*j represent, respectively, the deflection and slope of j. 

Equations 2 to 5 represent the boundary conditions whereas Equations 6 
and 7 are the initial conditions. Equation 7 simply states that at time 
zero, the member deflected shape is dependent on x, K, El and L. 

The first task toward the solution of the matrix equation is the 
assembly of the three coefficient matrices. The [K] matrix is assembled 
from the individual element matrices combined in such a way so as to enforce 
the given boundary and inter-element compatability conditions. To 
illustrate the procedure, an example of a beam with four elements as shown 
in Figure 2 is given in Appendix A. 

The global mass matrix is a diagonal form of a lumped mass matrix which 
was developed (Reference 6) for use with elements where translational 
degrees of freedom are mutually parallel, such as beam or plate elements. 
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This matrix may be written as: 

3 

2 

mL 2 

78 

m 

mL 2 

39 


[M] = 


where : 

m - mass at each degree of freedom - pL(A) 
p « mass density (mass/in 3 ) 

L - length of element (in) 

A =* cross sectional area of element (in 2 ) 

In order to calculate the damping matrix [C] , it is necessary to first 
determine the modal shape and natural frequencies of the system. This is 
accomplished numerically by solving the following eigen value problem using 
the Jacobi method (Reference 7): 

<[K] - w 2 [ M ] [ {$ } ] = {0} - (10) 

where : 

03 - natural frequency, 

{$ } - modal vector. 

Once 03 and { $} are known, determination of the damping matrix proceeds as 

/ 

described in Reference 7. 


m 

2 


mL 2 

78 


(9) 
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Once all three coefficient matrices have been assembled, the solution 
of Equation 1 may proceed using any one of several solution algorithms 
available . 

2 . 2 Newmark 1 s Method 

Newmark*s method for solving the dynamic equilibrium equation is 
sometimes called the trapezoidal method because it is based on a linear 
interpolation to find succeeding points. This is done by assuming: 


{D} t A - {D} t + A t{D) t 



{D} t + g{D) t+ A t 


and 


( 11 ) 


{i>) t+ A t - (D} t +a t ((i-y){D) t + yiD), 


( 12 ) 


where At is a time increment, and 3 and y are arbitrary constants. By 
substituting Equations 11 and 12 into Equation 1 written at time t = t + At, 
one gets (Reference 6) : 


JL 

3A^ 


[C] + 


• [M] (D} t+ 


At 




At 


+ 


[C] 


[M]i 


( {D) t + 1 (D) + (At) 3 - - 1 {D} t \ + 

l $At T 23 J 

/ / • / / \ N 

{D } t + — — { D } t + -A - lUD} t 


(13) 


3At z BAt \23 I 2 

For a known loading function we may solve Equation 13 for the deflection at 
time t - t + At using the deflection, velocity and acceleration at time t. 
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The algorithm for Newmark*s solution is as follows: 

1. Compute the coefficient matrices from geometric and material 
properties. 

2. At t = 0, set initial conditions by prescribing { D } t=0 and 

3. Use Equation 1 to solve for { D } t =0 . 

4. Solve Equation 13 for 

* 

5. Solve Equation 11 for {D} t+ ^ t . 

6. Solve Equation 12 for 

7. Set {D} t = {D} t+At ; (D} t - (D) t+At ; {D> t - {D} t+At . 

8. If t < total time desired, go to Step 4. 

9. Stop. 

This method of solution is unconditionally stable if Y>0.5 and 
$ > (2Y + 1) 2 /16. With y— 0.5 and 3 - 0.25, there are no amplitude errors 
in any sine wave motion regardless of its frequency, although the periods 
are overestimated. The mode shapes of the member in this study, however, 
are not known exactly. Nevertheless, Y and 3 values of 0.5 and 0.25 
respectively, were tentatively chosen. The suitability of these values is 
evaluated later in Section 4. 

The initial static deflection vector required in Step 2 of the 
algorithm may be determined using any one of the several classical 
structural analysis techniques. An approximate shape function for the 
member due to a specified midpoint displacement Ao at time t = 0 is taken in 
the following form (Reference 2) : 


A 1 


TTXj 

sin + 

L 


kL 

4ttEI 



(14) 
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where : 


“ o 

Aj - (15) 

kL 

1 + — 

2ttEI 


The initial slope of the member at any point is found by differentiating 
Equation 14 resulting in: 


J 


Tj TTXj k ZTTXj 

— cos + sin ■■ 

L L 2EI L 


(16) 


where Xj is the position of node j along the member length. 

2.3 Central Difference Formulation 

The governing equations and formulation of the coefficient matrices to 
be used in the central difference method of solution are precisely the same 
as those previously given for Newmark's method. Once these geometric and , 
physical properties are determined, one proceeds by writing the central 
difference expressions for both velocity and acceleration at an arbitrary 
time t: 


{D) t 


(D>t 


2 (At) 
1 

(At) 2 


< D >t +A t ’ U» t - 


A l 


< D >f A t * 2 < D >t + < D >t- 


At I 


(17) 


(18) 


Equations 17 and 18 may then be substituted into Equation 1 to yield, after 
s ome rearrangement : 



= < R >t 




[C] 

2 (At) 


(D} t . 


A 1 


( 19 ) 
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The initial conditions (D} 0 and { D } 0 are prescribed and (D} 0 is found by 
solving Equation 1. Once these are known Equations 17 and 18 may be solved 
simultaneously to yield the displacements (D)_At required to start the computations. 


(At) 2 

{D) t - (D) 0 - At{D> 0 + {D) 0 

U 2 


( 20 ) 


The solution algorithm for central difference is as follows: 

1. Compute the coefficient matrices from geometric and material 
properties . 

2. Set At « time step increment. 

3. Set initial conditions by prescribing {D} t=0 and { D > t =0 . 

4. Solve Equation 1 for {D) t=0 . 

5. Solve Equation 20 for { D } _^ t . 

6. Solve Equation 19 for {D) t+ ^ t . 

7. Set (D) t . At = { D ) t , and {D) t = {D) t+At . 

8. If t < total time desired, go to 6. 

9. Stop. 

The central difference method is a conditionally stable, explicit 
method of solution. Conditionally stable implies that if A t is not chosen 
small enough, the predicted response of the system will grow unbounded. A 
preliminary numerical study showed that At must be in the range from 0.001 
to 0.005, therefore, a ^t ® 0.001 sec. is used in this study. 
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3 . EXPERIMENTAL STUDY 


3.1 Specimen and Connection Details 

3.1.1 Specimen 

The experimental study consisted of conducting natural and forced 
vibration tests on a tubular aluminum member. The tests were performed both 
with and without passive damping devices present inside the member. The 
tubular member used was 14* -9” long with an outside diameter of 2” and a 
wall thickness of 0.125" , yielding an inside diameter of 1.75". A schematic 
of the member tested is shown in Figure 1. Note that the member was 
horizontal for all testing, with gravitational forces acting in the plane of 
motion. 

3.1.2 Connection Details 

The prototype end connection used in this study is shown in Figure 3. 

It is constructed of an aluminum alloy, weights 0.595 lbs. excluding 
fastener bolts, and has a volume of 3.988 in 3 . The connection has a total 
of nine clevis blades, six of which are in the horizontal plane. One of the 
blades is in the vertical plane (at C) and two are at 45 degrees to the 
horizontal plane. These two are located at 45 degrees relative to the 
vertical clevis and in the planes containing the two lower clevis blades 
shown in Figure 3(a). 

The fastner locations for the clevis blades in the horizontal plane are 
numbered 1 through 12. The member was fastened at locations 3 and 4 shown 
in Figure 3(a). Fasteners at locations 5 through 11 are used to mount the 
connection to a fixed base plate. No fastener was installed at location 12 
due to an interference problem with the support underneath the base plate. 
This did not make any difference since the other fasteners provided 
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sufficient fixity. Each fastener has a diameter of 0.25 in. and a length of 
0.94 in.. Washers were used at locations 1 and 2 only. 

The ends of the tubular member were threaded to allow one -half of the 
"snap- lock" connection to be screwed onto it. Small holes were drilled 
through this threaded connection and pins inserted to prevent rotation and 
loosening of the connection during testing. The other end of the snap -lock 
connection had its blade end fit snugly into one of the clevis blades of the 
prototype end connection and fastened by two bolts. The spring stiffness, 
k, shown in Figure 1 was determined by a statical analysis using an 
experimentally determined midspan deflection for a known concentrated load. 
This value was 53.1 k- in/rad. The assembled connection is shown in Figure 
4. 

3.2 Passive Damping Concepts 

Three different types of passive dampers referred to in Section 1.4 are 
described in this section. 

3.2.1 Copper Brush Dampers 

Figure 5 shows a copper brush damper 0.8125 inches in diameter, of 
total length 3.125 inches and a weight of 13.0 grams. The brush is 
manufactured by Oraack Industries, Onalaska, Wisconsin 54650 with a US 
Patent 41986 and an inventory control number 07668341989. It has a threaded 
aluminum piece 1.0 inch long at one end with a twisted wire 2.125 inches in 
length attached to it. The copper bristles are attached to the entire 
length of the twisted wire. This type of brush is commonly used in cleaning 
the bore of a 12 gauge shotgun. 

Figures 6 and 7 show schematically the attachments for the passive 
dampers and their spacing inside the tubular member. As shown in Figure 6, 
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the assembly consists of several parts. First, a helical spring with a 
stiffness of 0.44 lb/in. is attached to the inside of the connection through 
a hook on the snap -lock connector as shown in Figure 8. A nylon line is 
tied to the other end of the spring and also connected to the first copper 
brush damper. The nylon line (sportfisher monofilament line manufactured by 
K-Mart Corporation, Troy, Michigan 48084, 8013.9, No. EPM-40, inventory 
control number 04528201391) used in this investigation has a 40 lb. 
capacity. A series of nylon line and dampers are attached along the member 
length until the other end of the tubular member is reached. The end of the 
nylon line is passed through a hole in the snap- lock connector and stretched 
by an amount of 2.0 inches in the longitudinal direction to induce nominal 
tension in helical spring. It is then secured to the vertical clevis at the 
support. The stretched helical spring is shown in Figure 9. The resulting 
passive damping assembly is aligned with the longitudinal axis of the 
tubular member due to the small amount of axial tension. No axial 
compression of the member is induced by the passive damping assembly on the 
tubular member since both ends of the nylon line are connected to the rigid 
supports. Since the nylon line is flexible, a significant portion of the 
stretching is due to elongation of the line itself with the remainder of the 
stretching taking place in the spring. The dampers are installed 
equidistant ly between the ends of the member. 

As a part of the present study, the effect of both number of brushes 
and presence or absence of tension on the nylon line, on member damping was 
examined. 

In addition to baseline experiments on the specimen with no damping 
devices, a total of ten different conditions were examined. Tests with 1, 
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2, 3, 5, and 7 brushes were conducted both with and without tension in the 
line . 

3.2.2 Wool Swab Dampers 

Figure 10 shows a wool swab damper with a 1.0 inch diameter, a total 
length of 3 inches and a weight of 7.1 grams. The wool swab is manufactured 
by Omark Industries, Onalaska, Wisconsin 54650 with a US patent 415838 and 
an inventory control number 076683422187. It has a threaded piece at one 
end with a twisted wire attached to it to which the wool swab is attached. 
The aluminum piece is 0.75 inches long while the wool swab has a length of 
2.125 inches. This type of brush is commonly used for cleaning 12 gauge 
shotguns. The dampers are mounted inside the tubular member as shown in 
Figures 6 and 7. Tests were carried out using 1, 2, 3, 5, and 7 
equidistantly spaced wool swab dampers. 

3.2.3 Silly Putty in Chamber Dampers 

The final device examined was the "silly putty" in chamber damper shown 
in Figure 11. It consists of a sphere approximately 0.75 inches in diameter 
made from silly putty placed inside a hollow cylindrical chamber. Silly 
putty is a trade name for an elasto-plastic material commonly used as a 
children's toy. It is manufactured by Binney and Smith, Inc., Easton, PA 
18042, with an inventory control number of 07166208006. The chamber is made 
from a 1.0 in. long piece of a "Bristole Pipe" (PVC-1120, Schedule 40, ASTM- 
D-1785, nominal 1 inch pipe) having an original outer diameter of 1.058 in. 
and a wall thickness of 0.15 in. Since the damping effect was assumed to be 
provided by the silly putty, two steps were taken to reduce the mass of the 
damper thereby improving its efficiency. First, the inside diameter is 
increased by machining it to 0.914 in. resulting in a wall thickness of 0.07 
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in. Its weight is further reduced by drilling a total of seven 0.25 in. 
diameter holes around its periphery half-way from its ends. The silly putty 
is held inside the chamber by means of a plastic wrap ("Saran Wrap") 
stretched over the ends of the chamber and held in place with tape. The 
silly putty is then free to bounce around inside the chamber. The total 
weight of the damper including the silly putty, PCV chamber, and plastic 
wrap is 7.4 gms . The dampers are mounted inside the tubular member as shown 
in Figures 6 and 7. Tests were conducted with a nominal tension in the 
spring and with no tension in the spring using 1, 2, 3, 5 and 7 equidistant 
silly putty in chamber dampers. An additional test was performed with 11 
equidistant dampers and a nominal tension in the spring. 

3.3 Test Setup and Procedures 

The instrumentation used in the tests consisted of a proximity probe, 
harmonic vibration devices and a deflection- time plotter. This section 
summarizes the test setup and procedures followed for all the experiments 
included in this report. 

Figure 12 shows a schematic of the member natural vibration test setup. 
A weight, W - 7.9 lb. was suspended at the member midspan by means of a 
cord, causing a total midspan deflection of 5/32 in. To induce natural 
vibration, the cord was cut with a pair of scissors, thereby releasing the 
member. The time dependent deflection at member midspan is recorded by 
means of a proximity probe shown in Figure 13, and connected to a 
deflection-time plotter. 

Figure 14 shows the member forced vibration setup, a schematic of which 
is shown in Figure 15. Forced vibration of the specimen was obtained using 
a vibrator (Model 203-25-DC) with an oscillator (Model TPO-25). The 
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vibrator applies a forcing function of the type: 

■j 

F(t) - F 0 cos (fit) ( 2 i) 

in which F 0 - 4 lb. , t - time, and Q - frequency of the forcing function. 

The applied frequency may be controlled using the oscillator. 

The forcing function F(t) is transmitted from the vibrator to the 
tubular member through a fabricated vibrator connector as indicated in 
Figure 14. The details of this mechanical connector are shown in Figure 16. 
It consists of three main segments PQ, QR and RU interconnected at Q and R 
by means of pins. End P is connected to the vibrator. The end U is 
connected to the lower part of a metal hose clamp provided around the 
tubular member at midspan as indicated in Figure 14. The parts QR and RU 
can be disengaged at R by pulling out the pin RS instantaneously in the RS 
direction as indicated by the arrow at S. A string attached at S is used to 
pull out the pin. Once the pin is pulled, the arm QR drops freely and the 
beam is free to vibrate without constraints. Both joints Q and R are well 
lubricated to reduce friction. The vibrator connector in the engaged and 
the disengaged positions is shown in Figures 17(a) and 17(b), respectively. 

A record is made of the deflection- time response of the member once the 
forcing function, F(t) , is removed. 

3.4 Test Results and Discussion 

In this section, the results from the member natural and forced 
vibration tests are presented and discussed. 

3.4.1 Natural Vibration 

All passive damping concepts were tested with natural flexural member 
vibration caused by releasing a weight at midspan as explained in Section 
3.3. The initial midspan deflection, A 0 > due to the suspended weight is 
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0.1563 in. A summary of the test results for the tubular member with no 
dampers as well as with wool swab, copper brush, and silly putty in chamber 
dampers is given in Table 1. The number of dampers, the total weight of the 
damping assembly, the damping ratio and the damping efficiency index are 
listed for each passive damping assembly. The logarithmic decrement method, 
as described in Reference 8, was used with the experimentally obtained 
deflection versus time plots to obtain the damping ratio. 

The calculation of the damping ratio for the natural vibration tests 
was obtained using the first sixteen cycles and reading the amplitudes 
directly from the experimental deflection versus time plots. Each £ value in 
Table 1 was then obtained by taking the average results of three tests for 
each combination of damping devices. 

The efficiency index is defined (Reference 1 and 2): 


( 22 ) 


M. 


in which £ is the damping ratio with the damping devices, 4 0 is the damping 
ratio in the absence of any passive damping device, and M d is the total mass 
of the damping assembly. 

The natural frequency from all of the experiments was found to be 8.4 
Hz. The deflection versus time plots referenced in this section are 
obtained using the average £ value and natural frequency from the 
experiments, and the following A -t relationship (Reference 1). 


A 


sin a),,t + cos &) H t 


The damped circular frequency,U) d , is given by: 
ui/l - ~ 2 


(23) 


w d = w\Jl - (24) 

The details including the listing of a computer program utilizing Equation 
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23 to produce a deflection versus time plot are given in Reference 1. A 
baseline plot of deflection versus time for the member with no dampers is 
shown in Figure 18. 

3.4. 1.1 Copper Brush Dampers 

For the copper brush dampers the maximum £ - 0.0131 is obtained with an 
assembly of three damping devices. This assembly produces the maximum r) “ 
16.72 in/lb -sec 2 . Figure 19 shows the corresponding average A-t plot for a 
10 second duration. Figure 20 shows the effect of the three copper brushes 
on the deflection time envelopes. The vertical ordinate in this figure is 
designated by A e to indicate that the figure represents the envelopes rather 
than the complete A -t relationship. The damping ratios from the experiments 
are given in Table 2(a). In addition to the test conducted as described in 
Section 3.3, a series of tests were made with no tension in the damping 
assembly. These tests, conducted with 1, 2, 3, 5 and 7 devices in the 
specimen showed no significant increase in member damping regardless of the 
number of devices used. The results are summarized in Table 2(b). One 
plausible explanation for this is as follows. The outer diameter of the 
copper brush is less than the inside diameter of the member. When there is 
no tension in the damping assembly, the devices are free to bounce inside 
the specimen. Because the vibrations are relatively small and the natural 
frequency low, the assembly with no tension has a tendency to move with the 
specimen, bouncing slightly inside the member. Due to the relatively 
negligible mass of the damper as compared to the member this nearly 
coincident movement produces minimal damping of the vibration. With a 
slight tension in the assembly, it can have its own natural frequency, 
different from the specimen. As a result, when vibration of the specimen is 
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induced, the impact of the damping assembly with the side of the tube sets 
the assembly in motion. Two types of motion then contribute to the damping. 
First, because of the difference in natural frequency of vibration impact of 
the dampers against the inside of the tubular member acts to damp the 
vibration. Secondly, the frictional interaction between the dampers and the 
member inside surface takes place while the dampers vibrate both in plane 
but out of phase, and axially. When the number of dampers is increased 
beyond three with nominal tension, the damping ratio decreases. 

3 . 4 . 1 . 2 Wool Swab Dampers 

For the wool swab dampers the maximum £ = 0.0105 was obtained with an 
assembly of three dampers resulting in an efficiency of 9.05 in/lb-sec 2 . 

The maximum r) « 12.34 was obtained with a single damper assembly yielding a 
damping ratio of 0.0099. Figures 21 and 22 represent the A-t plots for the 
member with three, and one wool swab damper assemblies, respectively for a 
10 second duration. Figures 23 and 24 show the effects of these damping 
assemblies on the deflection- time envelopes. The damping ratio increased as 
the number of dampers was increased from one to three. Increasing the 
number of devices beyond three resulted in a decrease in both damping ratio 
and efficiency. The small negative efficiency noted for seven devices can 
be taken as practically zero. It was found that a variation in the method 
of attachment of the assembly to test specimen from concentric to an 
eccentric connection had no significant effect on the resulting damping 
ratio. The results are given in Tables 3(a) and 3(b). 

3.4. 1.3 Silly Putty in Chambers Dampers 

For silly putty in chambers dampers, the maximum 5 - 0.0115 was 
obtained with an assembly of three dampers resulting in a rj - 15.73 in/lb- 
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« 


sec 2 , whereas the maximum n *= 21.35 in/lb-sec 2 was obtained with an assembly 
of two dampers corresponding to Z, = 0.0113. Figures 25 and 26 represent the 
A-t plots for the member with three and two silly putty in chamber damper 
assemblies, respectively, for 10 second duration. Figures 27 and 28 show 
the effects of these damping assemblies on the deflection- time envelopes. 

The damping ratio was found to increase as the number of dampers was 
increased from one to three. Increasing the number of dampers beyond three 
resulted in a decrease of both damping ratio and efficiency. The tests 
conducted with no tension in the assembly showed a slight increase in 
damping ratio up to the three damper assembly. An increase in the number of 
dampers beyond three with no tension on the assembly showed no increase in 
damping ratio above the baseline damping ratio for the empty member. The 
results are given in Tables 4(a) and 4(b). Of all the passive damping 
devices tested in this study, the assembly of three silly putty in chamber 
dampers was found to be the most efficient. Therefore, these dampers were 
chosen for further study under forced harmonic vibration. 

3 . 5 Forced Then Free Vibration 

It was discovered during testing that the vibration employed for the 
forced vibration tests allowed only a limited amount of travel. This meant 
that the deflection of the member at the location where the vibrator was 
attached was limited to what the vibrator would allow. Nevertheless, forced 
vibration tests were conducted on the individual member since it was not 
known initially whether or not the dynamic deflection would exceed the 
vibrator capacity. The results presented later in this section indicated 
that the vibrator constrained the member deflection for a certain range of 
forcing function frequencies including that which would otherwise have 
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constituted a resonance condition. This limitation must be taken into 
account when evaluating the performance of the dampers on an individual 
member. 

3. 5. 1.1 Silly Putty in Chamber Dampers 

The results of the experimental study of the member under forced then 
free vibration are summarized in Table 5. Tests were conducted with no 
dampers, and 1, 2, and 3 dampers inside the member. Each of these 
assemblies was subjected to a force of 4 lb. at the member midspan, at 
frequencies of 2 , 5, 7, and 9 Hz, corresponding to ratios of 0.238, 

0.596, 0.834, and 1.073, respectively. An additional test was conducted on 
the empty member and the 3 damper assembly using a frequency of 12 Hz (ft/ 

« 1.430). The experimental results are shown in Figures 29 through 32. The 
free vibration part of the deflection- time graph is obtained by disengaging 
the forcing function from the member midspan as described in Section 3.3. 

The constrained dynamic deflection amplitude, A * D » and its dimensionless 
value, A*d/As> where As t * ie calculated static midspan deflection due to a 
4 lb. load, are listed in Table 5. The constrained dynamic deflection 
amplitude is the measured amplitude of the initial constrained force part of 
the deflection- time plots. Also listed in Table 5 are the maximum initial 
free vibration amplitudes, A F , for each assembly and frequency considered. 
Two dimensionless quantities are derived from this value as A F /A S and A F /A* D . 

The data in Table 5 shows that the A*d/As va lues range from 0.59 to 
0.95. For all the cases, the maximum value was observed for an applied 
force frequency of 5 Hz. It was also found that theA F /As an< * 

A f /A*d ratios were gradually increasing for increasing forcing function 
frequencies. One important consequence of the deflection constraint imposed 
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by the vibrator is that no resonance phenomenon could be produced in the 
vicinity of 8.4 Hz. The average damping ratios were obtained from the free 
vibration part of the deflection- time curves and are listed in Table 5. As 
seen from this data, the single silly putty in chamber damper configuration 
provides the maximum decrease in free vibration amplitude. Another 
important observation to be made is that the q values in Table 5 are 
significantly less than the corresponding values for the same damping 
assemblies given in Table 1. This is attributable to the dependence of the 
damping ratio on the initial velocity which is considerably greater for the 
results reported in Table 5 than for those in Table 1. 

3.5 Comparison of Damping Efficiencies 

In Section 3.4, the efficiency index based on Equation 22 was computed 
for each damping device. The average values of n and the associated damping 
assembly weight for natural vibration were presented in Table 1. Figure 33 
shows the curves between n and the weight of dampers used in the natural 
vibration tests for various damping concepts. The silly putty in chamber 
dampers provided the most efficienct damping of the member. It is worth 
noting that all. of the curves have ascending and descending portions which 
define the maximum attainable damping efficiency. In general, an increase 
in damping assembly weight beyond 50 grams results in a decline in 
efficiency. 

Figure 34 shows the relationships between the damping efficiency and 
the number of damping devices for natural vibration using all three 
concepts. These curves also show that, in general, an assembly of more than 
three damping devices results in a decline in efficiency. This may indicate 
that the first and second mode shapes are dominating the dynamic response. 
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By applying the dampers to locations in the vicinity of maximum deflection 
for these mode shapes, the maximum efficiency was realized. Any increase in 
the number of dampers beyond three adds mass to the system, and is 
associated with a decrease in damping. 

The average damping efficiency indices for the forced then free 
vibration tests for 1, 2, and 3 silly putty in chamber dampers are given in 
Table 5. The maximum efficiency was obtained using one silly putty in 
chamber damper and a forcing function frequency of 5 Hz. No correlation 
between the maximum efficiency and the initial vibration amplitude was 
observed. However, the maximum average damping ratio for each device was 
found to occur near the theoretical resonance of the member (between 7 and 9 
Hz) in spite of the inability of the apparatus to allow the resonance to 
occur. 
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4. NUMERICAL STUDY 


4.1 Natural Vibration 

4.1.1 Finite Element versus Experiment 

The formulation and solution algorithmn using Newmark's method for 
computing the dynamic response of a beam was given in Section 2. In this 
section, a comparison is made of the deflection versus time relations from 
this finite element analysis to those obtained experimentally. 

A preliminary study showed that for At - 0.0001 sec.,* the central 
difference formulation described in Section 2.3 gave precisely the same 
results as Newmark’s method. Since Newmark’s method provides accurate 
results even with larger time steps, it was used to produce Figure 35 
through 43. Figure 35 shows a comparison of the finite element and 
experimental A-t plots for the member with no dampers. The solid line is 
the finite element solution and the dashed line is the experimental curve 
using a frequency of 8.4 Hz and the average damping ratios from Table 1. 
Figure 36 shows a comparison of the finite element and experimental A -t 
plots for the member with 3 copper brush dampers. In both of these curves, 
it can be seen that the period of the vibration obtained using finite 
elements is exagerated by approximately 32%. However, the amplitudes of the 
vibration are accurate to within 5%. 

4.1.2 Finite Element versus Finite-Difference 

The A-t curves representing the finite-difference solution are obtained 
using the computer program developed in Reference 2. Figures 37 and 38 show 
the comparison of the finite element and the finite-difference solutions for 
the member with no dampers and three copper brush dampers, respectively. 

The data for these plots is obtained from Table 1. As indicated in these 
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figures, the difference in the period calculated by these two methods is 
approximately 26%. However, the amplitudes of the vibration from the two 
analyses are within 3% of each other. Figure 39 is a comparison of the 
finite-element and finite difference solutions for a simply supported beam 
(k t = k 2 - 0) . Similar correlation is also observed for a fixed end beam 
(k t « k 2 = oo ) . In the presence of end connections of intermediate fixity, 
the two analyses provide somewhat differing results. 

4.2 Forced then Free Vibration 

In this section, curves obtained from the finite element solution for 
various forcing function frequencies are given. Also, a comparison of the 
theoretical solution to experimental results is made for both the member 
with no dampers and the member with one silly putty in chamber damper at a 
forcing function frequency of 2 Hz. 

Figure 40 shows the response using Newmark's methqd for a beam with no 
dampers and subjected to a 4 lb. force at a frequency of 2 Hz. After 1 
second, the forcing function is removed and the beam is allowed to vibrate 
freely. Figure 41 shows the response of the same system with a forcing 
function frequency of 6 Hz. This frequency corresponds to a frequency ratio 
°f 0.95, where co fe - 6.3 Hz is the natural frequency of the beam from 
the finite element solution. Clearly, this represents a nearly resonant 
condition as expected. After 1 second, the forcing function is removed and 
the member is allowed to vibrate freely. 

Figures 42 shows the finite element and experimental curves for the 
member with no dampers and subjected to a 4.0 lb. force at a frequency of 2 
Hz. Although the forced vibration portions of the two curves at Q « 2 Hz 
are quite similar, the free vibration amplitudes differ significantly. The 
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reasons for this difference may be as follows. In the experiment, the 

forcing function was terminated by pulling the pin RS from the vibrator 

connector shown in Figure 16. During the tiny time interval in which the 
pin was being pulled out, the contact and frictional forces involved in 
disengaging the segment QR from RU were unintentionally transferred to the 
members thereby retarding its initial amplitude in the free vibration range. 
Consequently, the ensuing envelope of the experimental free vibration A-t 
curve is considerably narrower than the theoretical one. Similar effects 
are observed in Figure 43 which shows the finite element and experimental 
results when one silly putty in chamber damper is used. 

At larger ft values such as those of the order of 6 Hz, the fc-t 

relations from the finite element analysis do not match the experimental 

ones even in the forced vibration range. This is primarily due to the 

constraints imposed by the vibrator on the maximum member deflecting as 
explained earlier in Section 3.5. 

4.3 Finite Element Analysis for Forced Vibration 

As mentioned earlier, the vibrator used in the experimental study 
constrained the motion of the member in the presence of a forcing function. 
As a result, the actual effect of passing damping could not be observed for 

this condition. Therefore, a numerical study was conducted to examine the 

effect of passive damping in the presence of the forcing function. In this 
section, the theoretical results showing both the extent of damping which 
would occur during the forced vibration and the effect of the dampers on the 
deflection- time envelopes are presented and discussed. Figure 44 shows the 
theoretical dynamic magnification factor (DMF) , A D /A gl versus the frequency 
ratio for damping ratios of 0.0094, 0.0131 and 0.50. The first two 
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values of the damping ratios were obtained from the member tests with no 
dampers, and 3 copper brush dampers, respectively. As can be seen in this 
figure, the copper brush dampers do not change the DMF appreciably for non- 
resonance frequency ratios. However, the dampers reduce the DMF by 
approximately 7% at resonance. 

Figure 45 shows the deflection versus time relationship for the member 
with no dampers and with three copper brush dampers, with a forcing function 
frequency of 6.35 Hz — 1.0) for one second, and allowed to vibrate 

freely thereafter. These curves show that the passive damping results in a 
member amplitude reduction in the forced vibration range, however, its most 
beneficial effect occurs during the free vibration. After 3 seconds of free 
vibration, the amplitudes of the member with dampers are approximately 40% 
less than those corresponding to the empty member. 
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5. CONCLUSIONS AND FUTURE RESEARCH 

5.1 Conclusions 

The following conclusions are drawn from the research conducted herein: 

1. The silly putty in chamber concept provides the maximum passive 

\ 

damping efficiency under member natural vibration, as compared to 
the copper brush or the wool swab concepts. 

2. The copper brush concept provides the largest damping ratio of the 
system under natural vibration. 

3. Due to the limitation of the vibrator used, the effectiveness of 
the passive damping concepts could not be evaluated until the 
forcing function was disengaged. 

4. Frictional and contact forces acting on the member during 
disengagement from the vibration apparatus caused a reduction of 
the ensuing free vibration member amplitude. 

5. The theoretical results indicate that in the presence of a forcing 
function, the passive damping devices provide the most effective 
damping in the vicinity of the resonant frequency. 

6. The theoretical results show that passive dampers are 
considerably more effective under member natural vibration than 
during forced vibration. 

7. Under natural vibration, the finite element solution results in 
periods which are nearly 30 percent greater than the experimental 
ones. However, amplitudes are reasonably accurate. The accuracy 
of the results is improved when the member ends are pinned or 
fixed. 
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5.2 Future Research 


The most successful passive damping concepts identified herein should 
be examined using forced vibration equipment which would allow investigation 
of their effectiveness at or near resonance. Attempts should be made to 
identify a means of disengaging an applied force without adversely affecting 
the dynamic response of the member. These tests should be conducted on both 


individual members and structure sub-assemblies. 
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TABLES 



Table 1. Member natural vibration test results for copper 
brush, wood swab, and silly putty in chamber dampers. 


PASSIVE 

DAMPING 

CONCEPT 

NUMBER OF 
DAMPERS 

WEIGHT OF 
DAMPING 
ASSEMBLY 
(GM) 

AVERAGE 

DAMPING 

RATIO 

DAMPING 

EFFICIENCY 

INDEX 

(INAB. - SEC 2 ) 






No Dampers 


0.00 

0.0094 

0.00 

Copper 

i 

13.0 

0.0098 

5.39 

Brush 

2 

26.0 

0.0107 

8.76 

Dampers 

3 

39.0 

0.0131 

16.72 


5 

65.0 

0.0129 

9.44 


7 

91.0 

0.0097 

0.48 

Wool 

1 

7.10 

0.0099 

12.34 

Swab 

2 

14.20 

0.0102 

9.87 

Dampers 

3 

21.30 

0.0105 

9.05 


5 

35.50 

0.0101 

3.46 


7 

49.70 

0.0091 

-0.99 

Silly 

1 

7.8 

0.0100 

13.48 

Putty 

2 

15.6 

0.0113 

21.35 

in 

3 

23.4 

0.0115 

15.73 

Chamber 

5 

39.0 

1 

0.0109 

6.74 

Dampers 

7 

54.6 

0.0097 

0.96 


11 

85.8 

f i 

0.0094 

0.00 






















Table 2(a); Damping ratios from natural vibration tests with copper 

brushes and no cord tension. 


Number of Devices 


^ 2 

C 3 







i 

0.0097 

0.0095 

0.0095 

0.0096 

2 

0.0097 

0.0097 

0.0097 

0.0097 

3 

0.0098 

0.0096 

0.0096 

0.0097 

5 

0.0094 

0.0096 

0.0094 

0.0095 

7 | 

0 . 0095 

0.0095 

l 

0.0096 

1 

0.0095 


Table 2(b). Damping ratios from natural vibration tests with copper 

brushes and nominal cord tension. 


Number of Devices 

r 

'O 

1 

£ 

2 

C 

3 

C AVG 






i 

0.0097 

0.0102 

0.0096 

0.0098 

2 

0.0105 

0.0109 

0.0108 

0.0107 

3 

0.0133 

0.0128 

0.0131 

0.0131 

5 . 

I 

0.0129 

0.0131 

0.0128 

0.0129 

7 

0.0096 

0.0098 

0.0095 

0.0096 


34 





















Table 3(a); Damping ratios for natural vibration tests with wool 
brushes and concentric cord support. 


Number of Devices 

* 


C 3 

5 

AVG 

i 

0.0098 

0.0097 

0.0098 

0.0098 

2 

0.0100 

0.0098 

0.0101 

0.0100 

3 

0.0107 

1 

0.0108 

0.0101 

0.0105 

5 

0.0096 

0.0100 

0.0097 

0.0098 

7 

i 

0.0094 

i 

0.0094 

i 1 

0.0094 

| 

0.0094 


Table 3(b). Damping ratios for natural vibration tests with wool 

brushes and eccentric cord support. 


Number of Devices 

S 

l 

e 

2 

C 

3 

C 

AVG 






i 

0.0099 

0.0099 

0.0098 

0.0099 

2 

0.0103 

0.0101 

0.0101 

0.0102 

3 

0.0105 

0.0105 

0.0105 

0.0105 

5 

0.0099 

0.0102 

0.0102 

0.0101 

1 

7 

0.0093 

0.0090 

0.0090 

0.0091 
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Table 4(a). Damping ratios for natural vibration tests with silly 

putty and no cord tension. 


Number of Devices 

k 


S 3 

^ AVG 






i 

0.0096 

0.0097 

0.0096 

0.0096 

2 

0.0099 

0.0099 

0.0102 

0.0100 

3 

0.0101 

0.0101 

0.0101 

0.0101 

5 

0.0095 

0.0095 

0.0093 

0.0094 

7 

0.0094 

0.0092 

i 

0.0096 : 

0.0094 


Table 4(b). Damping ratios for natural vibration tests with silly 

putty and nominal cord tension. 


Number of Devices 

C i 

^2 

^3 

CMM 






i 

0.0104 

0.0096 

0.0100 

0.0100 

2 

0.0112 

0.0115 

0.0113 

0.0113 

3 

0.0112 

0.0115 

0.0117 

0.0115 

5 

0.0109 

0.0109 

0.0109 

0.0109 

7 

0.0096 

1 

0.0099 

0.0096 

0.0097 

11 

0.0094 

0.0094 

0.0094 

0.0094 
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Table 5. Member forced then free vibration test results for silly putty 

in chamber dampers . 


Passive 

Forcing 

Damping 

Function 

Concept 

Frequency 

(Hz) 



Va* d 


Average 
Damping 
Ratio 
? Avg. 



No 

Dampers 5 


1 Silly 
Putty 
in 

Chambei 

Damper 


2 Silly 
Putty 
in 

Chambei 

Damper 


3 Silly 
Putty 
in 

Chambei 

Damper 



60 

53 

83 


0.030 

0.057 

0.077 

0.090 

0.082 


0 . 

00 

0 . 

00 

O 

0 . 

00 

00 

0 . 

92 

0 . 

39 

0 . 

.76 

0 . 

.67 

1 . 

.05 

0 . 

00 

0 . 

.71 

0 

.97 

1 

.14 

1 

.03 


0.45 0.0043 

0.91 0.0072 

1.10 0.0073 

1.23 | 0.0069 

i 

1.75 I 0.0025 


0.48 

0.0049 

0.80 ! 

j 

0.0069 

0.80 

0.0076 

1.28 

0.0080 



0.47 

0.0057 

0.76 i 

| 

0.0055 

1.12 j 

j 

0.0064 

1.32 

0.0082 

1.75 

0.0036 


: 6.7 

; (-3.4) 
3.4 

i 

I 12.4 


10.5 

-12.7) 

(-6.7) 

9.7 

8.2 































FIGURES-?. 



Partial 

Rotational F(t) 




Figure 1. Schematic of tubular member 



Figure 2. Example of finite element model for the beam. 
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tb) Edge View between qq and bb 


Figure 3. Some end connection details 






DRIdNAD PAGE IS 
.OS POOR QUALITY 



Figure 4. Member end connection 


> 
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Figure 5. Copper brush damper 
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Figure 6. Schematic of attachments for passive 
damper inside tubular member 
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Figure 8 • Helical spring attachment at end e 
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Figure 9. 


Stretched helical spring at end e 
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Figure 11. Silly putty in chamber damper 
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Figure 12 • Schematic of member natural vibration setup 







Figure 13. Proximity probe 
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Figure 14. Member forced vibration setup 





Figure 15 . Schematic of member forced vibration setup 
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(a) 


Vibrator connector 
in engaged position 


(b) Disengaged vibrator 
connector 




Figure 17 . 


Vibrator connector in engaged and 
disengaged positions 
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Figure 18. 




ampe r s 


(3 
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Figure 20. Effect of 3 copper brush dampers on A“t envelope. 
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Figure 22. Ave 



No Dampers 
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Figure 23. Effect of 3 wool swab dampers on A-t envelope. 


No Dampers 



Figure 24. Effect of 1 wool swab damper on A-t envelope. 



Putty 





Silly Putty 
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Figure 26. Average A~t plot for member with 2 silly putty in chamber dampers. 



No Dampers 
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Figure 27. Effect of 3 silly putty in chamber dampers on 
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Figure 28. Effect of 2 silly putty in chamber 
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1.0 2.0 3.0 4.0- SiO- 6.0 7.0 


T (sec) 


(a) 2 Hz 


1.0 2.0 3.0 4.0 5.0 6.0 7.0 

T (sec) 

(b) 5 Hz 



T (sec) 


(c) 7 Hz 


1.0 2.0 3.0 4.0 5.0 6.0 7.0 

T (sec) 

(d) 9 Hz 


Figure 29. Experimental A-t plots for member "constrained" forced 
then free vibration with no dampers and frequencies of 2, 5, 7 and 9 Hz. 
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mm 


it::; 




|Htl»p««|M«ttVl|«i 

tstmw.«T ,!, i 


II 


1.0 2.0 3.0 4.0- 5 . 0 - 6.0 7.0 

T (sec) 

. (a) 2 Hz 


1.0 2.0 3.0 4.0 5.0 6.0 7.0 

T (sec) 

(b) 5 Hz 


rForced 


rForced 


|HI^SiH3Si!!l2i!!!l:0t!lIIfli.H(I!P^|niniR9RVVRVVVRi|9iBHHB| 


1.0 2.0 3.0 4.0* 5.0* 6.0 7.0 

T (sec) 


1.0 2.0 3.0 4.0 3.0 0.0 /.U 

T (sec) 


(c) 7 Hz (d) 9 Hz 

Figure 30. Experimental A-t plot for member "constrained" forced then free 
vibration with 2 silly putty in chamber damper and forcing function frequencies 
of 2, 5, 7 and 9 Hz. 
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Figure 32. Experimental A -t plots for member "constrained" forced then 
free vibration with 3 silly putty in chamber dampers and forcing function 
frequencies of 2, 5, 7, and 9 Hz. 
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Si I I y Putty 
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WEIGHT OF DAMPERS (gm) 

Figure 33. Damping efficiency index versus weight curves for natural vibration tests. 
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NUMBER OF DRMPERS 








10th cycle experiment 

/- 10th cycle finite element 



Figure 36. Finite element versus average experimental A~t plots for member with 3 copper brush dampers 




10th cycle finite element 



Figure 37. Finite element versus finite difference A-t plots for member with no dampers 






10th cycle finite element 
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Figure 38. Finite element versus finite difference A-t plots for member with 3 copper brush dampers 


finite element 
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Figure 39. Finite element versus finite difference A“t plots for simply supported beam with no dampers 
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Figure 40. Finite element A _t plot for member with no dampers* and subjected to a 4.0 lb. force at 2 Hz for 1 second 
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Figure 42. Theoretical and experimental 
force at 2 Hz for 1 second on member with 
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Figure 43. Theoretical and experimental forced then free A-t plots for a 4.0 lb. 
force at 2 Hz for 1 second on member with no dampers. 
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APPENDICES 



APPENDIX A 


EXAMPLE OF FOUR ELEMENT BEAM STIFFNESS MATRIX 


In this appendix the procedure used to assemble the beam stiffness 
matrix using a beam composed of four elements is presented. 

The typical element stiffness matrices for Elements b and c as shown in 


Figure 2 are 

given as 

(Reference 

6): 



12EI 

6EI 

-12EI 

6EI 



~L r 


~\ 7 ~ 



4 El 

-6EI 

2EI 

(K] b>c “ 


L 

L 2 


' 


12EI 

-6EI 


Symmetric 

L 3 

L 2 





4EI 





L 


Since only planar motion is considered, axial effects are negligible and, 
therefore, not included in the element stiffness matrix. 

Derivation of the stiffness matrix for Element a as shown in Figure 2 
is as follows. The flexibility matrix for the element is given by: 


[F] = [H] t [F] cl [H] + [F] m + [F] c2 (A. 2) 

in which [H] is the equilibrium matrix given by: 


10 0 


[H] = 


OIL 
0 0 1 


(A. 3) 


and [F] cl represents the flexibility of the connection at end one, [F] m is 
the flexibility of the element itself and [F] c2 is the flexibility of the 
connection at end two. These are defined as follows: 
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E F 3ci “ 


[F], 


1 ' 

0 


L 

4EIk 

L 3 

L 2 

3EI 

2EI 

L 2 

L 

2EI 

El 


[F] c2 - 


0 

0 


0 

L 


4EI 


therefore the flexibility matrix in full can be written as 




[F] 


The inverse of [F] is given by: 



( K 2 2 ). 


[FJ 


-1 


-3EI (2k + 1) 


[ K u] a = 

[H] [K 22 ] [H] 1 


L 3 ” ' 

(k + 1) 

~l7~ 

(k + 1) 

-3EI 

(2k + 1) 

El 

(4k + 3) 


(k + 1) 

L 

(k + 1) 

matrices now follow 

from 

: 

3E1 

(4k + 1) 

3EI 

(2k) 

L 3- 

(k + 1) 

\7~ 

(k + 1) 

3EI 

(2k) 

El _ 

(4k) 


(k + 1) 


L (k + 1) 


86 _ 


(A. 4) 


(A. 5) 


( A . 6 ) 


(A. 7) 


(A. 8) 


(A. 9) 



c 


■ 87 



3EI (4k + 1) 
L 3 (1 + k) 

[K 12 ] d = [K 21 ] fc d = 

- [H] [K] 

-3EI (2k + 1) 
L 2 (1 + k) 


6EI k 
L 2 (1 + k) 


2EI (k) 

L (k + 1) 


(A. 14) 


[Kid 



(4k + 1) 3EI 
(1 + k) L r ~ 


El 

L 


Symmetric 


(2k + 1) 
(1 + k) 

(4k + 3) 
(1 + k) 


-3EI (4k + 1) 


L 3 

(k + 1) 

-3EI 

(2k + 1) 

L 2 

(1 + k) 

3EI 

(4 + 1) 

L 3 ” 

(1 + k) 


6EI 

(k) 

~L 2 ~ 

(k + 1) 

2EI 

(k) 

L 

(1 + k) 


-6EI 

(k) 

L 2 

(1 + k) 

4EI 

(k) 

L 

(1 + k) 


(A. 15) 


Using the above element matrices, the following global matrix is assembled: 



[*lJa 

[*22 1 a 

[0] 

[0] , 

[0] 




[*21 la 

[*2 2 ]a + [*ll] 

b [*l 2 lb 

[0] 

[0] 



[K] = 

[0] 

[*21 lb 

[*22l b +[*lll 

c [*l 2 lc 

[0] 


(A. 16) 


[0] 

[0] 

[*21 1c 

[*22lc + [*ll]d 

[*12 1 d 




[0] 

[0] 

[0] 

[*21 Id 

[*22 1 d 



This 

is an n 

x n matrix 

where n - 2N + 

2 , and N = the 

number 

of 

elements 


The first two boundary conditions are enforced by putting 1.0 in the 
diagonal corresponding to the translational degrees of freedom at the 
supports and setting all other entries in that row and column equal to zero. 
The last two boundary conditions are accounted for in the derivation of the 
individual stiffness matrices. 

Note that an adjustment to the stiffness matrix must be made when the 
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APPENDIX B 


COMPUTER PROGRAMS 

As a part of this study, two computer programs were developed to solve 
the dynamic equilibrium matrix equation given in Chapter 2. A brief 
description of these programs along with their listings and sample outputs 
are given in this appendix. 

B . 1 NEWMARK 

This program is based on the analysis described in Section 2.2. A 
description of the required input data is given at the beginning of the 
program listing. Data is input by means of the data statements in lines 48 
to 53 of the program listing. The output consists of the time in seconds 
and corresponding midspan deflection in inches. 

B . 2 CENDIF 

Program CENDIF is based on the analysis described in Section 2.3. Data 
input and output are the same on NEWMARK. 
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FILE: NEWMARK FORTRAN A OLD DOMINION UNIVERSITY — CMS — 4.2 


NEWOOOIO 

IMPLICIT REAL *8 (A-H.O-Z) NEW00020 

NEW00030 

DOUBLE PRECISION L,K1 ,K2,K(70,70) ,M(70,70) , C (70, 70) , U (70) ,UDT(70) ,NEW00040 
6 RT (70) ,CI (70,70) ,C2 (70,70) ,C9 (70,70) ,C4(70) ,KEL(4,4) ,F (70) ,C5(70) .NEW00050 
6 DUM (70) , UT (70) ,MINV(70,7O) ,C 8 (70,70) ,UDTP (70) ,UDDT(70) ,UDDTP (70) .UNEW00060 


STP (70) , UTN (70) , C3 (70) , C 6 (70) , FREQ (5) , K I NV (70, 70) , X (70, 70) , A (70,70) NEW00070 
6 , B (70,70) , D (70) , E I GV (70) , DAMRAT (70) ,ANS (70,70) , RES (70,70) ,XT (70, 7 ONEWOOO 8 O 



&) 

1 FPR=0 


NEW00090 

NEW00100 

C ************* * ft* * * ft* INPUT DATA ************* ** * ************ * 

NEW001 10 

c** 


** 

NEW001 20 

c** 

IFPR = PRINT REQUEST VARIABLE FOR JACOBI 

** 

NEW00130 

c** 

0 = 00 NOT PRINT INTERMEDIATE VALUES 

** 

NEW00140 

c** 

1 = PRINT INTERMEDIATE VALUES 

** 

NEW00150 

c** 


** 

NEW00160 

c** 

L = LENGTH (IN) 

** 

NEW00170 

c** 


** 

NEW00180 

c** 

NUMEL = NUMBER OF ELEMENTS (MUST BE AN EVEN NUMBER) 

** 

NEW00190 

c** 


** 

NEW00200 

c** 

TS = TIME STEP; DELTA 1 T 1 (SEC) 

** 

NEW002 1 0 

c** 


** 

NEW00220 

c** 

ROW = MASS PER UNIT LENGTH (K 1 P*SEC**2/ 1 N**2) 

** 

NEW00230 

c** 


** 

NEW00240 

c** 

E = MODULUS OF ELASTICITY (KSI) 

** 

NEW00250 

C** 


ft* 

NEW00260 

c** 

XI = MOMENT OF INERTIA (IN**4) 

ft* 

NEW00270 

c** 


ft* 

NEW00280 

c** 

AR = AREA (IN**2) 

ft* 

NEW00290 

c** 


ft* 

NEW00300 

c** 

K 1 = ROTATIONAL STIFFNESS AT END 1 (K* 1 N/RAD) 

ft* 

NEW00310 

c** 


ft* 

NEW00320 

c** 

K2 = ROTATIONAL STIFFNESS AT END 2 (K* 1 N/RAD) 

ft* 

NEW00330 

c** 


ft* 

NEW00340 

c** 

ZETA = DAMPING RATIO 

ft* 

NEW00350 

C** 


** 

NEW00360 

c** 

TT = TOTAL TIME FOR PROGRAM EXECUTION (SEC) 

*ft 

NEW00370 

c** 


ft* 

NEWOO 38 O 

C** 

PO = MAGNITUDE OF THE FORCING FUNCTION (KIPS) 

ft* 

NEW00390 

c** 


** 

NEW00400 

c** 

OMEGA = FREQUENCY OF THE FORCING FUNCTION (RAD) 

ft* 

NEW00410 

c** 


** 

NEW00420 

c** 

DELO = PRESCRIBED INITIAL DEFLECTION AT MEMBER MIDSPAN 

ft* 

NEW00430 

c** 

(IN) 

ft* 

NEW00440 

c** 


** 

NEW00450 

c * * ***** ***** ***** ********** ** * ** ****** ****** ** ** ** * * * ** ** * ****** 

DATA L,NUMEL,TS,ROW/ 177, 10,0.000500,202. 1454E-09/ 

DATA E,XI ,AR/ 10000., 0.32500000,0. 7363/ 

DATA K1,K2,TT/53.1100,53.1100,6.00/ 

DATA PO, OMEGA, ZETA/O. 002226367, 12.56637062,0.007200/ 

DATA GAMA, BETA, DELO/0. 50, 0.25, 0.0079 1 1 */ 

1 C0UNT=0 

NEW00460 

NEW00470 

NEW00480 

NEW00490 

NEW00500 

NEW00510 

NEW00520 

NEW00530 

NEW00540 

NEW00550 
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FILE: NEWMARK FORTRAN A OLD DOMINION UNIVERSITY — CMS — 4.2 



T 1 ME=0 .0 

NEW00560 

NEW00570 


WRITE (2,1) 

NEW00580 


H=L/NUMEL 

NEW00590 

1 

FORMAT (/IX, 'THIS IS NEWMARKS SOLUTION 1 ) 

NEW00600 


WRITE (2,1059) K 1 

NEW00610 

1059 

FORMAT (/IX, 'STIFFNESS -- ' , D 1 6 . 9) 

NEW00620 


WRITE (2, 1060) ZETA 

NEW00630 

1060 

FORMAT (/IX, 'DAMPING ' , D 1 6 . 9) 

NEW00640 


WRITE (2, 1061) PO 

NEW00650 

1061 

FORMAT (/IX, 'FORCE ' , D 1 6 . 9) 

NEWOO 66 O 


WRITE (2,1062) OMEGA 

NEW00670 

1062 

FORMAT (/IX, 'FREQUENCY — ' , D 1 6 . 9) 

NEWOO 68 O 


N=2*NUMEL+2 

NEW00690 

NEW00700 

C 

WRITE (2,176) 

NEW00710 

NEW00720 

C 

WRITE (2,177) 

NEW00730 


DO 10 1=1 , N 

NEW00740 

NEW00750 


DO 10 J=1 , N 

NEW00760 

10 

K ( 1 , J) =0 .0 

NEW00770 


K (1,0-1000 

NEWOO 78 O 

NEW00790 

NEW00800 


K (N- 1 , N- 1 ) =1000 

NEW00810 


K (2 , 2) =E*X 1 *4 . *K 1 / (H* (K 1 + 1 .)) 

NEW00820 

NEWOO 83 O 

NEW00840 


K (2 , 3) = (- 1 .)*3.*E*XI*2.*K1/((H**2)*(K1 + 1 .)) 

NEW00850 


K(2,4)=E*XI*2.*K1/(H*(K1+1 • )) 

NEWOO 86 O 


K (3,2) -K (2,3) 

NEW00870 


K (3 , 3) =3 . *E*X 1 * (4 . *K 1+1 .) / ( (H**3) * (K 1 + 1 . ) ) 

NEWOO 88 O 


K (3,4) = (-3.) *E*XI*(2.*K1 + 1 .) / ( (H**2) * (K1 + 1 .)) 

NEWOO 89 O 


K (4,2) =K (2,4) 

NEW00900 


K (4, 3) *K (3,4) 

NEW00910 


K (4, 4) =E*X 1 * (4 .*K 1+3-) / (H* (Kl + 1 .)) 

NEW00920 


K (N-3,N~3)=3.*E*XI*(4.*K2+1 .) /((H**3) *(K2+1 .)) 

NEW00930 

NEW00940 

NEW00950 


K (N-3, N-2) = (3 .) *E*X 1 * (2 . *K2+1 .) / ( (H**2) * (K2+1 .) ) 

NEW00960 


K(N-3,N)=6.*E*XI*K2/((H**2)*(K2+1 •)) 

NEW00970 


K (N-2,N~3)=K (N-3, N-2) 

NEW00980 


K (N-2 , N-2) =E*X 1 * (4 . *K2+3 .) / (H* (K2+ 1 . ) ) . 

NEW00990 


K (N-2 , N) =2 .*E*X 1 *K2/ (H* (K2+1 .)) 

NEW01000 


K (N,N-3)=K (N-3, N) 

NEW01010 


K (N , N-2) =K (N-2.N) 

NEW01020 


K (N , N) =E*X 1 *4 . *K2/ (H* (K2+1 .) ) 

NEW01030 



NEW01040 


KEL(l.l) = 12.*E*XI/(H**3) 

NEW01050 

NEW01060 

NEW01070 


KEL (1 ,2) = 6.*E*XI/(H**2) 

NEW01080 


KEL (1 , 3) = (-1 .) *KEL (1,1) 

NEW01090 


KEL (1 ,4) = KEL (1 , 2) 

NEW01 100 
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ILE: 

NEWMARK FORTRAN A OLD DOMINION UNIVERSITY -- CMS -- 4.2 

' 


KEL (2 , 2) = 4.*E*XI/(H) 

NEW01 1 10 


KEL (2 , 3) = (-l.)*KEL(l,2) 

NEW01 120 


KEL (2 , 4) = KEL (2 , 2) /2 . 

NEW01 130 


KEL (3 . 3) = KEL (1,1) 

NEW01 140 


KEL (3.4) -■ KEL (2.3) 

NEW01 150 


KEL (4,4) = KEL (2 , 2) 

NEW01 160 
NEW01 170 


IF (K (2,2) .LE. 0.00001) K (2 , 2) = 1 .0 

NEW01 180 


IF (K(N,N) .LE. 0.00001) K(N,N) = 1.0 

NEW01 190 


DO 30 1=1,4 

NEW01200 


DO 30 J=1 ,4 

NEW01210 

30 

IF (J . GT. 1 ) KEL (J , 1 ) =KEL (1 , J) 

NEW01220 

NEW01230 


DO 50 JK =2 , NUMEL- 1 

NEW01240 


1 1 =JK*2-2 

NEW01250 


JJ=I 1 

NEW01260 


DO 45 1=1,4 

NEW01270 


DO 40 J=1 ,4 

NEW01280 

NEW01290 


K ( 1 l + l , J J+J) = K (1 l + l , JJ+J) +KEL (1 , J) 

NEW01300 

NEW01310 

40 

CONTINUE 

NEW01320 

^5 

CONTINUE 

NEW01330 

50 

CONTINUE 

NEW01340 * 
NEW01350 

66 

DO 75 1-1, N 

NEW01360 


DO 70 J = 1 , N 

NEW01370 


M ( 1 , J) =0.0 

NEW01380 


C (1 , J)=0.0 

NEW01390 

70 

CONTINUE 

NEW01400 

75 

CONTINUE 

NEW01410 


c 

NEW01420 


M (1 , 1) =39- 

NEW01430 


M (2 , 2) =H**2 

NEW01440 


M (N- 1 , N- 1 ) =39 • 

NEW01450 


M(N,N)=H**2 

NEW01460 



NEW01470 


DO 80 l=3,N-3,2 

NEW01480 


J=l + 1 

NEW01490 


M ( 1 , 1 ) =78 . 

NEW01500 


M(J,J)=2.*(H**2) 

NEW01510 

80 

CONTINUE 

NEW01520 

NEW01530 


DO 90 1=1, N 

NEW01540 

90 

M ( 1 , l)=M(l , l)*(ROW*H/78.) 

NEW01550 


CALL JACOBI (K,M,N,IFPR,X,EIGV) 

NEW01560 


DO 95 1-1, N 

NEW01570 

95 

DAMRAT ( 1 ) =ZETA 

NEW01580 

NEW01590 


CALL DAMP(N,EIGV,X,M, DAMRAT, C) 

NEW01600 

NEW01610 

NEW01620 


PRINT*, 1 IN START' 

NEWOI63O 



NEW01640 


DO 300 1=1, N 

NEW01650 
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300 

RT ( 1 ) =0 . 0 

NEW01 660 


PI=ACOS (-1.0) 

NEW01670 


RT (N/2) =P0* (DCOS (OMEGA*T 1 ME) ) 

NEW01680 

NEW01690 


CALL INVERT (M,MINV,N) 

NEW01700 

NEW01710 


DO 333 l-l.NUMEL+l 

NEW01720 


DUM1=K1*L/ (4*PI*E*XI) 

NEW01730 


. Z»PI*H*(I-1)/L 

NEW01740 


UT (2* 1 ) = (DELO/ (1 .+DUM1 *2) ) * ( (P 1 /L*DCOS (Z) ) + (DUM1 *2*P 1 /L*DS 1 N (2*Z) 

NEW01750 


&)) 

•NEW01760 


UT (2* 1 - 1 ) = (DELO/ (1 .+DUM1 *2) ) * (DS 1 N (Z) +DUM1 * (1 . -DCOS (2*Z) ) ) 

NEW01770 

333 

CONTINUE 

NEW01780 


DO 302 1=1, N 

NEW01790 


SUM=0 .0 

NEW01800 


DO 301 J=1 , N 

NEW01810 

301 

SUM=SUM+K ( 1 , J) *UT (J) * (- 1 . ) 

NEW01820 

302 

DUM ( 1 ) =SUM 

NEW01830 
NEW01 840 


DO 306 1=1, N 

NEW01850 


SUM=0 .0 

NEW01860 


DO 305 J=1 ,N 

NEW01870 

305 

SUM=SUM+M 1 NV ( 1 , J) *DUM (J) 

NEW01880 

306 

UDDT ( 1 ) =SUM 

NEW01890 

NEW01900 


DO 310 1=1 ,N 

NEW01910 


UDT ( 1 ) =0 .0 

NEW01920 

310 

CONTINUE 

NEW01930 

NEW01940 


PRINT*, 'OUT START' 

NEW01950 

NEW01960 


DO 320 1=1, N 

NEW01970 


DO 320 J=1,N 

NEWOI 98 O 

320 

C 1 ( 1 , J) = K ( 1 , J) +GAMA*C ( 1 , J) / (BETA*TS) +M ( 1 , J) / (BETA* (TS**2) ) 

NEW01990 


CALL INVERT (Cl ,C2,N) 

NEW02000 

120 

DO 122 1=1, N 

NEW02010 


C3 (1 ) = (GAMA/ (BETA*TS) ) *UT ( 1 ) + (GAMA/BETA- 1 .0) *UDT ( 1 ) +TS* ( (GAMA/ (BE 

NEW02020 


&TA*2 .) ) -1 .0) *UDDT ( 1 ) 

NEW02030 

NEW02040 


C4 ( 1 ) =UT ( 1 ) / (BETA* (TS**2) ) +UDT ( 1 ) / (BETA*TS) + ( ( 1 . / (2 . *BETA) ) - 1 . 0) *UNEW02050 


6DDT ( 1 ) 

NEW02060 

122 

CONTINUE 

NEW02070 

NEW02080 


DO 130 1=1, N 

NEW02090 


SUM=0 .0 

NEW02100 


DO 125 J=1 ,N 

NEW021 10 

125 

SUM=SUM+C(I ,J)*C3(J) 

NEW02120 

130 

C5(0=SUM 

NEW02130 



NEW02140 


DO 410 1=1, N 

NEW02150 


SUM=0 .0 

NEW02160 


DO 400 J=1,N 

NEW02170 

400 

SUM=SUM+M ( 1 , J) *C4 (J) 

NEW02180 

410 

C6 (1 ) =SUM 

NEW02190 


NEW02200 
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DO 4 20 1=1 ,N NEW02210 

420 F ( I ) =RT ( I ) +C5 ( I ) +C 6 ( I ) NEW02220 

NEW02230 

DO 430 1=1 , N NEW02240 

SUM=0.0 NEW02250 

DO 425 J=1 * N NEW02260 

425 SUM=SUM+C2 (I ,J)*F (J) NEW02270 

430 UTP ( I ) =SUM NEW02280 

NEW02290 

I COUNT= I COUNT+1 NEW02300 

TIME=TIME+TS NEW02310 

IF (TIME. GT. 1 .00) GO TO 199 NEW02320 

RT (N/2) =P0* (DCOS (OMEGA*T I ME) ) NEW02330 

GO TO 200 NEW02340 

199 RT (N/2) =0.0 NEW02350 

200 EXACT=0 .0 NEW02360 

JEST=1 NEW02370 

I F ( I COUNT. EQ. 10) GO TO 1 4 1 . NEW02380 

GO TO 143 NEW02390 

141 WRITE (2, 175) TIME, UTP (N/ 2 ) , JEST NEW02400 

I C0UNT=0 NEW02410 

NEW02420 

143 DO 150 1=1, N NEW02430 

UDDTP ( I ) = (UTP ( I ) -UT ( I ) - (TS*UDT ( I ) ) - ( (TS**2) * (0 . 5 -BETA) *UDDT ( I ) ) ) *NEW02440 
6 ( 1 . / ( (TS**2) *BETA) ) NEW02450 

150 CONTINUE NEW02460 

NEW02470 

DO 161 1=1, N NEW02480 

UDTP ( I ) =UDT ( I ) +TS* ( ( ( 1 . 0-GAMA) *UDDT ( I ) ) + (GAMA*UDDTP ( I ) ) ) NEW02490 

UT ( I ) =UTP ( I ) NEW02500 

UDT ( I ) =UDTP ( I ) NEW02510 

UDDT (I ) =UDDTP (I ) NEW02520 

161 CONTINUE NEW02530 

IF (TIME. GT .TT) GO TO 500 NEW02540 

NEW02550 

GO TO 120 NEW02560 

175 FORMAT (F10.8, 1X.F10.8, IX, II) NEW02570 

176 FORMAT (/IX , 1 TIME DEFLECTION AT L/2 ' ) NEW02580 

177 FORMAT (' NEW02590 

£ • ) NEW02600 

500 STOP NEW02610 

END NEW02620 

NEW02630 

NEW02640 

SUBROUTINE I NVERT (AO , A , N) NEW02650 

DOUBLE PRECISION A (70 , 70) , AO (70 , 70) NEW02660 

NEW02670 

DO 1 1=1, N NEW02680 

DO 1 J=1 , N NEW02690 

1 A ( I , J) =A0 ( I , J) NEW02700 

NEW02710 

NP=N+1 NEW02720 

A (1 ,NP) =1.0 NEW02730 

DO 10 1=2, N NEW02740 

10 A (I ,NP) =0.0 NEW02750 
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DO 40 J=1 ,N 

NEW02760 

NEW02770 


DO 20 LX=1 , N 

NEW02780 

20 

A (NP , LX) =A (1 ,LX+1)/A (1 , 1) 

NEW02790 


DO 30 KX=2 , N 

NEW02800 


DO 30 LX=1 , N 

NEW02810 

30 

A (KX-1 ,LX) =A (KX.LX+l) -A (KX, 1) *A (NP.LX) 

NEW02820 


DO 40 LX=1 , N 

NEW02830 

40 

A (N , LX) =A (NP.LX) 

NEW02840 


RETURN 

NEW02850 

NEW02860 


END 

NEW02870 


SUBROUTINE JACOBI (K ,M, N , 1 FPR, X , E 1 GV) 

NEW02880 

C 

SUBROUTINE JACOBI 

NEW02890 


IMPLICIT REAL*8 (A-H.O-Z) 

NEW02900 


DOUBLE PRECISION A (70 , 70) , B (70 , 70) , X (70 , 70) , E 1 GV (70) , D (70) , 

NEW02910 

&K (70,70) ,M (70, 70) 

NEW02920 


1 FPR=0 

NEW02930 

C 

COMMON/K.M/ 

NEW02940 

C 

WRITE (2, 1051) 

NEW02950 

C1051 

FORMAT (/IX, ' INPUT DATA ') 

NEW02960 

C 

READ (1, *) N, 1 FPR 

NEW02970 

C 

WRITE (2 , 1001) N, 1 FPR 

NEW02980 

C 

DO 1010 1=1, N 

NEW02990 

C 

READ (1 ,*) (A (1 ,J) ,J=1 ,N) 

NEW03000 

C 

WRITE (2,1110) (A ( 1 , J) , J=1 , N) 

NEW03010 

C1010 

CONTINUE 

NEW03020 

C 

DO 1020 1=1, N 

NEW03030 

C 

READ (1 ,*) (B (1 ,J) ,J=1 , N) 

NEW03040 

C 

WRITE (2,1110) (B ( 1 , J) ,J=1,N) 

NEW03050 

C1020 

CONTINUE 

NEW03060 

C1001 

FORMAT (2 110) 

NEW03070 

Cl 1 10 

FORMAT (8F10.4) 

NEW03080 


DO 2 1=1, N 

NEW03090 


DO 1 J= 1 , N 

NEW03100 


A ( 1 , J) =K ( 1 , J) 

NEW031 10 


B (1 , J) =M ( 1 , J) 

NEW03120 

1 

CONTINUE 

NEW03130 

2 

CONTINUE 

NEW03140 


NSMAX=15 

NEW03150 

C 

WRITE (2,1980) 

NEW03160 

1980 

FORMAT (/IX, 1 EIGENVALUES ') 

NEW03170 


RT0L=1 . D- 1 2 

NEW03180 


1 0UT=2 

NEW03190 


DO 10 1=1, N 

NEW03200 


1 F (A (1 , 1) .GT.O.AND.B (1 ,1) .GT.O.) GO TO 4 

NEW03210 


WRITE (1 OUT, 2020) 

NEW03220 


STOP 

NEW03230 

4 

D(I)=A(I, l)/B(l,l) 

NEW03240 

10 

E 1 GV ( 1 ) =D (1 ) 

NEW03250 


DO 30 1=1 ,N 

NEW03260 


DO 20 J=1 , N 

NEW03270 

20 

X (1 , J)=0. 

NEW03280 

30 

X ( 1 , 1 ) = 1 . 0 

NEW03290 


IF (N . EQ. 1 ) RETURN 

NEW03300 
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C 



NEW03310 

C IN 

ITIALIZE SWEEP COUNTER AND EIGEN ITERATION 


NEW03320 

C 



NEW03330 


NSWEEP =0 


NEW03340 


NR=N-1 


NEW03350 

40 

NSWEEP=NSWEEP+1 


NEW03360 


IF (IFPR.EQ.l) WRITE (1 OUT, 2000) NSWEEP 


NEW03370 


PRINT *, 1 SWEEP NUMBER... ' .NSWEEP 


NEW03380 

C 



NEW03390 

C CHECK IF PRESENT OFF DIAGONAL ELEMENT IS TOO 

LARGE 

NEW03400 

C 



NEW03410 


EPS= (0.01 **NSWEEP) **2 


NEW03420 


DO 210 J= 1 ,NR 


NEW03430 


JJ=J+1 


NEW03440 


DO 210 K 1 =JJ , N 


NEW03450 


1 F (DABS (A (J , Kl) ) . LT . 1 . D-20) GO TO 211 ‘ 


NEW03460 


EPTOLA= (A ( J , K 1 ) *A ( J , K J ) ) / (A (J , J) *A (K 1 , K 1 ) ) 


NEW03470 


GO TO 212 


NEW03480 

211 

EPT0LA=0 .0 


NEW03490 

212 

EPTOLB= (B ( J , K 1 ) *B ( J , K 1 ) ) / (B (J , J) *B (K 1 , K 1 ) ) 


NEW03500 


IF ((EPTOLA.LT.EPS) .AND. (EPTOLB.LT. EPS))GO TO 210 

NEW03510 


AKK=A (Kl , K 1) *B (J , Kl ) -B (K 1 , K 1) *A (J , Kl ) 


NEW03520 


AJJ=A(J,J)*B(J,K1) -B(J,J)*A (J,K1) 


NEW03530 


AB=A (J , J) *B (K 1 , K 1 ) -A (K 1 , K 1 ) *B (J , J) 


NEW03540 


CHECK= (AB*AB+4 . *AKK*AJ J) /4 . 


NEW03550 


IF (CHECK) 50,60,60 


NEW03560 

50 

WRITE (1 OUT, 2020) 


NEW03570 


STOP 


NEW03580 

60 

SQCH=DSQRT (CHECK) 


NEW03590 


D 1 =AB/2 .+SQCH 


NEW03600 


D2=AB/2 . -SQCH 


NEW03610 


DEN=D1 


NEW03620 


IF (DABS (D2) .GT.DABS (D 1 ) ) DEN=D2 


NEW03630 


IF (DEN) 80,70,80 


NEW03640 

70 

CA=0. 


NEW03650 


CG=(-1 . ) *A (J , Kl ) /A (Kl , K 1 ) 


NEW03660 


GO TO 90 


NEW03670 

80 

CA=AKK/DEN 


NEW03680 


CG=(-1 .) *AJ J/DEN 


NEW03690 

90 

IF (N-2) 100, 190, 100 


NEW03700 

100 

JP1=J+1 


NEW03710 


JM1=J~ 1 


NEW03720 


KP 1 =K 1+1 


NEW03730 


KM1=K 1 - 1 


NEW03740 


IF (JM1-1) 130, 110, 1 10 


NEW03750 

1 10 

DO 120 1 = 1 , JM1 


NEW03760 


AJ=A (1 ,J) 


NEW03770 


BJ=B (1 ,J) 


NEW03780 


AK=A (1 ,K1) 


NEW03790 


BK=B (1 , K 1 ) 


NEWO 38 OO 


A ( 1 , J) =AJ+CG*AK 


NEWO 38 IO 


B (1 , J) =BJ+CG*BK 


NEW03820 


A (1 , Kl) =AK+CA*AJ 


NEW03830 

120 

B (1 ,K1)=BK+CA*BJ 97 


NEW03840 

130 

IF (KP1-N) 140,140,160 


NEW03850 
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1 40 DO 150 l-KPl.N 
AJ=A (J, I) 

BJ=B (J, I) 

AK=A (K 1 , I ) 

BK-B(Kl.l) 

A (J , I ) =AJ+CG*AK 
B (J , I ) =BJ+CG*BK 
A (K 1 , l)=AK+CA*AJ 
150 B (K 1,1) =BK+CA*BJ 
160 IF (JP1-KM1) 170,170,190 
170 DO 180 I-JP1..KM1 
AJ=A (J , I ) 

BJ=B (J, I) 

AK=A (I , K 1 ) 

BK=B (I , K 1 ) 

A (J, I ) =AJ+CG*AK 
B (J , I ) =BJ+CG*BK 
A (I , Kl) =AK+CA*AJ 
180 B (I ,K1) =BK+CA*BJ 
190 AK=A (K 1 , K 1) 

BK=B (Kl , Kl) 

A (K 1 , Kl ) =AK+2 . *CA*A (J , K 1 ) +CA*CA*A (J , J) 

B (Kl ,K1)=BK+2.*CA*B (J , K 1) +CA*CA*B (J , J) 

A (J , J) =A (J, J)+2.*CG*A (J , K 1 ) +CG*CG*AK 
B (J , J) =B (J , J) +2 . *CG*B (J , K 1) +CG*CG*BK 
A (J ,K1) = 0 . 

B (J,K1) = 0 . 

C 

C UPDATE EIGENVECTOR MATRIX 
C 

DO 200 1-1, N 
XJ=X ( I , J) 

XK-X(I.Kl) 

X (I , J)=XJ+CG*XK 
200 X ( I ,K1)=XK+CA*XJ 
210 CONTINUE 
C 

C UPDATE EIGENVALUES 
C 

DO 220 1-1, N 

IF (A (I , I) . GT.O.AND.B (I , I) . GT.O) GO TO 220 
WRITE (IOUT, 2020) 

STOP 

220 EIGV (I)— A (I , l)/B(l , I) 

IF (I FPR.EQ.O) GO TO 230 
WRITE (I OUT, 2030) 

WRITE (IOUT, 2010) (E I GV ( I ) , I = 1 , N) 

CHECK FOR CONVERGENCE 

230 DO 240 1-1, N 
TOL=RTOL*D(l) 

D I F=DABS (EIGV(I)-D (I)) 

IF (DIF.GT.TOL)GO TO 280 
240 CONTINUE 


NEWO 386 O 
NEW03870 
NEWO 388 O 
NEWO 389 O 
NEW03900 
NEW03910 
NEW03920 
NEW03930 
NEW03940 
NEW03950 
NEWO 396 O 
NEW03970 
NEW03980 
NEW03990 
NEW04000 
NEW04010 
NEW04020 
NEW04030 
NEW04040 
NEW04050 
NEW04060 
NEW04070 
NEW04080 
NEW04090 
NEW04100 
NEW041 10 
NEW041 20 
NEW04130 
NEW04140 
NEW04150 
NEW04 1 60 
NEW04170 
NEW04180 
NEW04190 
NEW04200 
NEW04210 
NEW04220 
NEW04230 
NEW04240 
NEW04250 
NEW04260 
NEW04270 
NEW04280 
NEW04290 
NEW04300 
NEW043I0 
NEW04320 
NEW04330 
NEW04340 
NEW04350 
NEW04360 
NEW04370 
NEW04380 
NEW04390 
NEW04400 
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C 

C 

C 


CHECK ALL OFF D I AG ELEMENTS TO SEE IF ANOTHER SWEEP IS REQ'D 

EPS=RT0L**2 
DO 250 J=1 ,NR 
JJ=J+1 

DO 250 K 1 =J J , N 

IF (DABS (A (J , Kl) ) .LT.1.D-30)GO TO 251 
EPSA= (A (J , K 1 ) *A (J , K 1 ) ) / (A (J , J) *A (K 1 , K 1 ) ) 

GO TO 252 

251 EPSA=0.0 

252 EPSB= (B (J , Kl) *B (J , K 1) ) / (B (J , J) *B (Kl , K 1) ) 

IF ((EPSA.LT.EPS) .AND. (EPSB.LT.EPS)) GO TO 250 
GO TO 280 
250 CONTINUE 


C 

C 

C 


FILL OUT BOTTOM TRIANGLE OF RESULTANT MATRICES 6 SCALE EIGENVECTORS 


255 


260 


270 

C 

C 

C300 

310 


DO 260 1=1, N 
DO 260 J=1 ,N 
A (J , I ) =A ( I , J) 

B (J , I ) =B ( I ,J) 

DO 270 J=1,N 
BB=DSQRT (B (J , J) ) 

DO 270 Kl=l ,N 
X (K 1 , J) =X (K 1 , J) /BB 
WRITE (IOUT, 310) 

DO 300 1=1, N 

WRITE (I OUT, 20 10) (X (I ,J) ,J=1 ,N) 

FORMAT (/IX, ' THE EIGENVECTORS ARE ') 


C 

C 

C 


RETURN 


UPDATE THE ’ D ' MATRIX AND START NEW SWEEP IF ALLOWED 


280 

290 


2000 

2010 

2020 


DO 290 1=1 ,N 
D (l)-EIGV(l) 

IF (NSWEEP.LT. NSMAX) GO TO 40 
GO TO 255 

FORMAT (/IX, 1 SWEEP NUMBER IN JACOBI = 1 , 1 4) 
FORMAT (/1X.6E20.12) 


FORMAT (/IX, 
6DEFINITE ') 
2030 FORMAT (/IX, 
END 


:’c * 5’c ERROR SOLUTION STOP / MATRICES NOT POSITIVE 

CURRENT EIGENVALUES IN JACOBI ARE ') 


NEW044 10 
NEW04420 
NEW04430 
NEW04440 
NEW04450 
NEW04460 
NEW04470 
NEW04480 
NEW04490 
NEW04500 
NEW04510 
NEW04520 
NEW04530 
NEW04540 
NEW04550 
NEW04560 
NEW04570 
NEW04580 
NEW04590 
NEW04600 
NEW046 10 
NEW04620 
NEW04630 
NEW04640 
NEW04650 
NEW04660 
NEW0467O 
NEW04680 
NEW04690 
NEW04700 
NEW04710 
NEW04720 
NEW04730 
NEW04740 
NEW04750 
NEW04760 
NEW04770 
NEW04780 
NEW04790 
NEW04800 
NEW04810 
NEW04820 
NEW04830 
NEW04840 
NEW04850 
NEW04860 
NEW04870 
NEW04880 
NEW04890 
NEW04900 
NEW04910 


SUBROUTINE DAMP (N , E I GV, X,M, DAMRAT, C) 

IMPLICIT REAL*8(A-H,0-Z) 

DOUBLE PRECISION X(70,70) ,T(70,70) ,M(70,70) ,C (70,70) ,EIGV(70) , DAMRNEW04920 
£AT (70) NEW04930 

NEW04940 

DO 10 1=1, N NEW04950 
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E l GV ( I ) =DSQRT (E I GV ( I ) ) 

DO 10 J=1 ,N 
10 C(l,J)=0.0 

DO 20 I 1 = 1 ,N 

DA=2.*DAMRAT (I l)*EIGV(l I) 

DO 20 1=1, N 
DO 20 J=1 , N 

20 C ( I , J) -C ( I , J) +X ( I , I I ) *X (J , I I ) *DA 

DO 30 1=1, N 
DO 30 J=1 ,N 
T ( I ,J)=0.0 
DO 30 K1=1,N 

30 T ( I , J) =T ( I , J) +M ( I , K 1 ) *C (K 1 , J) 

DO 40 1=1, N 
DO 40 J=1 , N 
C (I , J) =0.0 
DO 40 K 1 = 1 ,N 

40 C(l ,J)=C(I.,J)+T(I ,K1)*M(K1,J) 

C DO 50 1=1, N 

C50 WRITE (2, 120) (C (I ,J) ,J=1 ,N) 

120 FORMAT (6D 14. 4) 

RETURN 

END 


NEW04960 
NEW04970 
NEW04980 
NEW04990 
NEW05000 
NEW05010 
NEW05020 
NEW05030 
NEW05040 
NEW05050 
NEW05060 
NEW05070 
NEW05080 
NEW05090 
NEW05100 
NEW051 10 
NEW05120 
NEW05130 
NEW05140 
NEW05150 
NEW05160 
NEW05170 
NEW05180 
NEW05190 
NEW05200 
NEW05210 
NEW05220 
NEW05230 
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OUT 


A OLD DOMINION UNIVERSITY — CMSL 4.0 8706 


THIS IS NEWMARKS SOLUTION 


STIFFNESS — 0.53H00006D+02 

DAMPING 0.130999982D-01 

FORCE 0.400000066D-02 

FREQUENCY — 0 . 251 300049D+02 

TIME DEFLECTION AT L/2 


0.00500000 0.07560742 1 
0.01000000 0.07411487 1 
0.01500000 0.07822479 1 
0.02000000 0.07755833 1 
0.02500000 0.07321689 1 
0.03000000 0.07439081 1 
0.03500001 0.07596175 1 
0.04000001 0.07030083 1 
0.04500001 0.06480695 1 
0.05000001 0.06403831 1 
0.05500001 O.O5838368 1 
0.06000001 0.04628549 1 
0.06500001 0.03770042 1 
0.07000001 0.03044851 1 
0.07500001 0.01549482 1 
0.08000001 -0.00155219 1 
0.08500001 -0.01381828 1 
0.09000001 -0.02866951 1 
0.09500002 -0.04917751 1 
0.10000002 -0.06667677 1 
0.10500002 -0.08058453 1 
0.1 1000002 -0.09807491 1 

0.1 1500002 -0.1 1601054 1 
0.12000002 -0.12788046 1 
0.12500002 -0.13797468 1 
0.13000002 -0.14964225 1 
0.13500002 -0.15657816 1 
0.14000002 -0.15730855 1 
0.14500002 -0.15770900 1 
0.15000002 -0.15621028 1 
0.15500003 -0.14753221 1 
0.16000003 -0.13531903 1 
0.16500003 -0.12327482 1 
0.17000003 -0.10672159 1 
0.17500003 -0.08461688 1 
0.18000003 -0.06275240 1 
0.18500003 -0.04085973 1 
0.19000003 -0.01462197 1 
0.19500003 0.01269540 1 
0.20000003 0.03673648 1 
0.20500003 0.06085743 1 
0.21000003 0.08652612 1 
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IMPLICIT REAL *8 (A-H.O-Z) 


CEN00010 

CEN00020 

CEN00030 

c 

THIS PROGRAM SOLVES FOR THE DEFLECTIONS OF A BEAM 


CEN00040 

CEN00050 

c 

SUBJECT TO A FORCING FUNCTION AT THE MIDPOINT 


CEN00060 

CEN00070 

c 

USING CENTRAL DIFFERENCE METHOD 

• 


CEN00080 
CEN00090 
CEN00100 
CEN001 10 


DOUBLE PRECISION L , K 1 , K 2 , K (70 , 70) , M (70 , 70) , C (70 , 70) 

,U (70) , UTN (70) 

, CENOO 1 20 


6 RT ( 70 ) , C 1 ( 70 , 70 ) ,C 2 ( 70 , 70 ) ,C 3 ( 70 , 70 ) , C4 ( 70 ) ,KEL (4,4) ,B( 70 ) ,C5(70) .CEN 00130 


6 DUM (70) , UT (70) , Ml NV (70,70) ,C 6 (70,70) ,X (70,70) ,EIGV(70) , DAMRAT (70) 
6 UDDT ( 70 ) 

1 fpr=o 

, CENOO 1 40 
CEN00150 
CEN00160 

C * * * * * * * * 5*c 5*c * * * * * * * * * * s ET CONSTANT DATA *********************** 

CEN00170 

c** 


** 

CEN00180 

c** 

L = LENGTH (IN) 

** 

CEN00190 

c** 


** 

CEN00200 

c** 

NUMEL = NUMBER OF ELEMENTS (MUST BE AN EVEN NUMBER) 

** 

CEN002 10 

c** 


** 

CEN00220 

c** 

TS * TIME STEP; DELTA ' T ' (SEC) 

** 

CEN00230 

c** 


** 

CEN00240 

c** 

ROW = MASS PER UNIT LENGTH (KIP*SEC**2/IN**2) 

** 

CEN00250 

c** 


** 

CEN00260 

c** 

E = MODULUS OF ELASTICITY (KSI) 

** 

CEN00270 

c** 


** 

CEN00280 

c** 

XI = MOMENT OF INERTIA (IN**4) 

** 

CEN00290 

c** 


** 

CEN00300 

c** 

AR « AREA ( 1 N**2) 

** 

CEN00310 

c** 


** 

CEN00320 

c** 

K 1 = ROTATIONAL STIFFNESS AT END 1 (K*l N/RAD) 

** 

CEN00330 

c** 


** 

CEN00340 

c** 

K2 = ROTATIONAL STIFFNESS AT END 2 (K* 1 N/RAD) 

** 

CEN00350 

c** 


** 

CEN00360 

c** 

ZETA = DAMPING RATIO 

** 

CEN00370 

c** 


** 

CENOO 38 O 

c** 

TT ■ TOTAL TIME FOR PROGRAM EXECUTION (SEC) 

** 

CEN00390 

c** 


** 

CEN00400 

c** 

PO = MAGNITUDE OF THE FORCING FUNCTION (KIPS) 

** 

CEN00410 

c** 


** 

CEN00420 

c** 

OMEGA = FREQUENCY OF THE FORCING FUNCTION (HZ) 

** 

CEN00430 

c** 


** 

CEN00440 

c** 


** 

CEN00450 

c**************************************************************** 

DATA L.NUMEL.TS, ROW/1 77.. 12,0.00010, 181 .95270-09/ 

DATA E,XI .AR/10000., .325, .7363/ 

DATA K 1 , K 2 ,TT/ 00 . 000 , 00 . 000 , 4 .00/ 

DATA PO, OMEGA, ZETA/O. 00, 0.000000000, 0.0000/ 

1 CQUNT=0 

CEN00460 

CEN00470 

CEN00480 

CEN00490 

CEN00500 

CEN00510 

CEN00520 

CEN00530 

CEN00540 

CEN00550 
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TIME=0.0 

CEN00560 

CEN00570 


H=L/NUMEL 

CEN00580 



CEN00590 


N=2*NUMEL-2 

CEN00600 

CEN00610 


DO 10 1=1, N 

CEN00620 


DO 10 J=1 , N 

CEN00630 

10 

K (1 , J) =0.0 

CEN00640 

CEN00650 


K (1 , 1) =3 • *E*X 1 * (4.*K1+1 .)/( (H**3) * (K 1 + 1 .)) 

CEN00660 


K ( 1 , 2) = (- 1 . ) *3 . *E*X 1 * (2 . *K 1+1 .) / ( (H**2) * (K 1+1 . ) ) 

CEN00670 


K (2, 1) =K (1 ,2) 

CENOO 68 O 


K (2,2)=E*XI*(4.*Kl+3.) /(H*(K! + 1 .) ) 

CEN0O69O 



CEN00700 


K (N- 1 ,N-1) =3.*E*X 1 * (4 . *K2+1 .) / ( (H**3) * (K2+1 .) ) 

CEN00710 


K (N- 1 , N) =3 . *E*X 1 * (2 . *K2+1 .) / ( (H**2) * (K2+1 . ) ) 

CEN00720 


K (N.N-l) =K (N-l ,N) 

CEN00730 


K (N , N) =E*X 1 * (4 . *K2+3 •) / (H* (K2+1 .)) 

CEN00740 

CEN00750 


KEL (1 , 1) = 12.*E*XI/(H**3) 

CEN00760 


KEL (1 ,2) = 6 .*E*X 1 / (H**2) 

CEN00770 


KEL (1,3) = (-1 .) *KEL (1,1) 

CEN00780 


KEL (1,4) = KEL (1,2) 

CEN00790 


KEL (2 , 2) = 4 .*E*X 1 / (H) 

CEN00800 


KEL (2,3) = (“1 .) *KEL (1 ,2) 

CEN00810 


KEL (2,4) = KEL (2,2) /2. 

CEN00820 


KEL (3 » 3) = KEL (1,1) 

CEN00830 


KEL (3.4) = KEL (2 , 3) 

CEN00840 


KEL (4,4) = KEL (2,2) 

CEN00850 

CEN00860 


DO 30 1=1,4 

CEN00870 


DO 30 J=1 ,4 

CENOO 88 O 

30 

IF (J . GT . 1 ) KEL (J , 1 ) =KEL ( 1 , J) 

CEN00890 

CEN00900 


DO 50 JK =1 , NUMEL-2 

CEN00910 


ll=JK*2-2 

CEN00920 


JJ=I 1 

CEN00930 


DO 45 1=1,4 

CEN00940 


DO 40 J— 1 ,l» 

CEN00950 

CEN00960 


K (1 l + l , JJ+J) = K (1 l + l ,JJ+J)+KEL (1 , J) 

CEN00970 

CENOO 98 O 

40 

CONTINUE 

CEN00990 

45 

CONTINUE 

CEN01000 

50 

CONTINUE 

CEN01010 

CEN01020 

C*******yoV***********^***********^*/c*fc***yc**^c*y<**3Sc**********A********A**CEN0 1 030 

CEN01040 

C 

THIS IS A TEST OF THE STIFFNESS MATRIX FOR THE STATIC LOAD CASE 

CEN01050 

CEN01060 

C 

DO 61 1=1, N 

CEN01070 

C6l 

RT ( 1 ) =0 . 0 

CEN01080 

C 

RT (N/2) = . 10 

CEN01090 

C 

CALL INVERT (K,C6,N) 

CEN01 100 
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c 

DO 63 1-1, N 


CEN011 10 

c 

SUM=0 .0 


CEN01120 

c 

DO 62 J=1 , N 


CEN01 1 30 

C 62 

SUM=SUM+C6 ( 1 , J) *RT (J) 


CEN01 140 

C 63 

U (l)=SUM 


CEN01 150 

c 

DO 65 1-1, N 


CEN01 160 

c 

PRINT*, 'DEFLECTION AT NODE ',1,' IS 

1 ,U(I) 

CEN01 170 

C 65 

WRITE (2,64) 1 ,U(I) 


CEN01 180 

C64 

FORMAT (/IX, 'STATIC SOLUTION... U(',l 

2,') = ',023.16) 

CEN01 190 

C 

GO TO 500 


CEN01200 

C 



CENO 1210 

C )'c * * * * * Vc ;V * Vc * * * * * * i'c :'c ft * * * * * * ft ft ft ft ft * * * ft * * ft ft ft ft * * ft ft * * ft * ft * * ft ft * ft ft ft * * * * ft ft * * 

CENO 1220 




CEN01230 

66 

DO 75 1-1, N 


CEN01240 


DO 70 J-l.N 


CEN01250 


M ( 1 ,J)=0.0 


CEN01260 


C (1 ,J)=0.0 


CEN01270 

70 

CONTINUE 


CENO 1280 

75 

CONTINUE 


CENO1290 

CEN01300 

CEN01310 


DO 80 1 =1 , N- 1 , 2 


CEN01320 


J=l + 1 


CEN01330 


M ( 1 , 1 ) =78 . 


CEN01340 


M(J, J)=2.*(H**2) 


CENO1350 

80 

CONTINUE 


CEN01360 

CEN01370 


DO 90 1-1, N 


CEN01380 

90 

M ( 1 , 1 ) =M ( 1 , 1) * (ROW*H/78.) 


CEN01390 

CEN01400 

C 

DO 100 1 = 1, N 


CENO 1410 


1 FPR=0 


CEN01420 


CALL JACOBI (K,M,N, 1 FPR.X.E IGV) 


CEN01 430 
CENO 1 440 


DO 100 1=1, N 


CEN01 450 

100 

DAMRAT (1 ) =ZETA 


CEN01 460 
CEN01470 


CALL DAMP (N,EIGV,X,M, DAMRAT, C) 


CEN01480 
CEN01490 
CEN01500 
CENO 1510 

C100 

C (1 , 1 ) =ZETA 


CEN01520 

CEN01530 

C***** PRINT STIFFNESS, MASS, AND DAMPING 

MATRICES *********** 

CEN01540 

C 



CEN01550 

c 

WRITE (2,220) N/2 


CEN01560 

c 

DO 210 1-1, N 


CEN01570 

C210 

WRITE (2,215) (K (1 , J) , J-l , N/2) 


CEN01580 

C 

WRITE (2 , 22 1 ) N/2 


CEN01590 

C 

DO 211 1=1, N 


CEN01600 

C211 

WRITE (2,215) (K (1 , J) ,J=N/2+l,N) 


CEN01610 

C 

WRITE (2,222) N/2 


CEN01620 

c 

DO 212 1-1, N 


CEN01630 

C212 

WRITE (2,215) (M (1 , J) ,J=1 ,N/2) 


CEN01640 

C 

WRITE (2,223) N/2 


CEN01650 
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C 

DO 213 1=1, N 

CEN01660 

C2 1 3 

WRITE (2,215) (M(l ,J) , J=N/2+l , N) 

CEN01670 

C 

WRITE (2 ,22k) N/2 

CEN01680 

C 

DO 2 1 k 1=1, N 

CEN01690 

C214 

WRITE (2,215) (C(l ,J) ,J=1 , N/2) 

CEN01700 

C 

WRITE (2 ,225) N/2 

CEN01710 

c 

00 216 1=1, N 

CEN01720 

C2 1 6 

WRITE (2,215) (C ( 1 ,J) ,J=N/2+l,N) 

CEN01730 

CEN01740 

215 

FORMAT (/IX, 5D 14.7) 

CEN01750 

CEN01760 

220 

FORMAT (/IX, 'THESE ARE THE FIRST ',13,' COLUMNS OF K ') 

CEN01770 

221 

FORMAT (/IX, 'THESE ARE THE LAST ',13,' COLUMNS OF K ') 

CEN01780 

222 

FORMAT (/IX, 'THESE ARE THE FIRST ',13,' COLUMNS OF M ') 

CEN01790 

223 

FORMAT (/IX, 'THESE ARE THE LAST M3,' COLUMNS OF M ') 

CEN01800 

22k 

FORMAT (/IX, 'THESE ARE THE FIRST ',13,' COLUMNS OF C ') 

CEN01810 

225 

FORMAT (/IX, 'THESE ARE THE LAST ',13,' COLUMNS OF C ') 

CEN01820 

CEN01830 

CEN018U0 


PRINT*, ' IN START' 

CEN01850 

CEN01860 


DO 300 1=1, N 

CEN01870 

300 

RT (1) =0.0 

CEN01880 

CEN01890 

C 

RT (N/2) =P0* (DS 1 N (OMEGA*T 1 ME) ) 

CEN01900 

CEN01910 


CALL INVERT (M,MINV,N) 

CEN01920 

CEN01930 

CEN01940 


P 1 =A COS (-1 .0) 

CEN01950 


DO 333 1=1 ,NUMEL 

CEN01960 


z=pi*h*i/l 

CEN01970 


UT (2*1) = (.1563/(1 .+DUM1*2) ) * ( (PI/L*DCOS (Z) ) + (DUM1 *2*P 1 /L*OS 1 N (2*Z) CEN01 98O 


&) ) 

CEN01990 


UT (2*1 -1) = (. 1563/(1 .+DUM1 *2) ) * (DSI N (Z)+DUM1* (1 .-DCOS (2*Z))) 

CEN02000 

333 

CONTINUE 

CEN02010 

C 

DO 334 1=1. N 

PRINT*, 'UT(' , 1 , ') = ' ,UT(I) 

CEN02020 

C334 

CEN02030 

CEN02040 



CEN02050 

CEN02060 


DO 302 1=1, N 

CEN02070 


SUM=0.0 

CEN02080 


DO 301 J-l.N 

CEN02090 

301 

SUM=SUM+K ( 1 , J) *UT (J) * (-1 .0) 

CEN02100 

302 

DUM ( 1 ) =SUM 

CEN021 10 
CEN02120 


DO 306 1=1, N 

CEN02130 


SUM=0.0 

CEN02140 


DO 305 J-l.N 

CEN02150 

305 

SUM=SUM+M 1 NV ( 1 , J) *DUM (J) 

CEN02160 

306 

UDDT ( 1 ) =SUM 

CEN02170 

CEN02180 


DO 303 1=1, N 

CEN02190 

303 

UTN (1) =UT (1) +UDDT ( 1 ) * (TS**2) /2 . 

CEN02200 
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CEN02210 

CEN02220 

c 

DO 302 1 = 1, N 

CEN02230 

c 

SUM=0.0 

CEN02240 

c 

DO 301 J=1 ,N 

CEN02250 

C301 

SUM=SUM+MINV(I ,J)*RT(J) 

CEN02260 

C302 

UTN ( 1 ) =SUM 

CEN02270 

CEN02280 

C 

DO 303 1=1, N 

CEN02290 

C303 

UTN ( 1 ) =UTN ( 1 ) * (TS**2) / (2 . 0) 

CEN02300 

CEN02310 


DO 310 1=1, N 

CEN02320 


U(l)=0.0 

CEN02330 

c 

UT ( 1 ) =0 . 0 

CEN02340 

310 

CONTINUE 

CEN02350 

CEN02360 


PRINT*,' OUT START' 

CEN02370 

CEN02380 


DO 320 1=1, N 

CEN02390 


DO 315 J-l.N 

CEN02400 


C 1 ( 1 , J) = M ( 1 , J) / (TS**2) + C ( 1 ,J)/(2.0*TS) 

CEN02410 


C2 (1 , J) = K ( 1 , J) - (2 . *M ( 1 , J) / (TS**2) ) 

CEN02420 


C3(I,J) = (M ( 1 , J) / (TS**2) ) - (C (1 ,J)/(2.*TS)) 

CEN02430 

315 

CONTINUE 

CEN02440 

320 

CONTINUE 

CEN02450 

CEN02460 


CALL INVERT(C1,C6,N) 

CEN02470 

CEN02480 

c 

DO 122 1=1, N 

CEN02490 

C 122 

WRITE (2,123) 1 ,C6 (1 , 1) 

CEN02500 

123 

F0RMAT(/1X, 'THIS IS C6 ( ' , 1 2 , ' ) ' , D23 • 1 6 ) 

CEN02510 

C 

WRITE (2, 176) 

CEN02520 

C 

WRITE (2, 177) 

CEN02530 

130 

DO 440 1=1, N 

CEN02540 


SUM=0 .0 

CEN02550 


DO 410 J=1 , N 

CEN02560 

C 

WRITE (2, 13DTIME.C2 (1 , J) , UT ( 1 ) 

CEN02570 

MO 

SUM=SUM+C2 ( 1 , J) *UT (J) 

CEN02580 

440 

C4 (1) =SUM 

CEN02590 

131 

FORMAT (/IX, 'TIME ' , F5 - 3, ' C2 ' , Dl8. 10, ' UT ',018.10) ' 

CEN02600 


DO 442 1=1, N 

CEN02610 


SUM=0.0 

CEN02620 


DO 411 J=1,N 

CEN02630 

41 1 

SUM=SUM+C3 (1 ,J) *UTN (J) 

CEN02640 

442 

C5(I)=SUM 

CEN02650 

CEN02660 


RT (N/2) =P0* (DS 1 N (OMEGA*T 1 ME) ) 

CEN02670 

CEN02680 


DO 140 1=1, N 

CEN02690 

140 

B(I)=RT(I)-C4(I)-C5(I) 

CEN02700 

CEN02710 

C 

DO 142 1=1, N 

CEN02720 

Cl 42 

WRITE (2, 149) 1 ,B(I) 

CEN02730 

149 

FORMAT (/IX, 'THIS IS B ( ' , 1 2 , ' ) = ',D23.l6) 

CEN02740 

CEN02750 
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C WRITE (2, 176) CEN02760 

C WRITE (2, 177) CEN02770 

DO 542 1=1 ,N CEN02780 

SUM=0.0 CEN02790 

DO 51 1 J= 1 ,N CEN02800 

511 SUM=SUM+C 6 (I ,J)*B(J) CEN02810 

542 U ( I ) =SUM CEN02820 

CEN02830 

I COUNT=l COUNT+1 CEN02840 

T I ME=T I ME+TS CEN02850 

CEN02860 

C SUM=0 .0 CEN02870 

C DO 199 1=1, 7, 2 CEN02880 

C X=l CEN02890 

C D 1 = ( (E*X I / (ROW* (L**4) ) ) **0.5) * (X**2) * (9-869604404) CEN02900 

C D2= ( (DS I N ( (X) *3. 14592654/2.0) )**2) CEN02910 

C D3=l .0- (DCOS (D 1 *T I ME) ) CEN02920 

Cl 99 SUM=SUM+D2*D3/(D1**2) CEN02930 

C199 SUM=SUM+D3/(D1**2) CEN02940 

C EXACT= (2 . *P0/ (ROW*L) ) *SUM CEN02950 

CEN02960 

C D 1=2 . *P0* (DS I N (OMEGA*TI ME) ) * (L**3) / ( (3 . 1 4592654**4) *E*X I ) CEN02970 

C SUM=0.0 CEN02980 

C DO 199 1=1 ,NUMEL-1 ,2 CEN02990 

C X= I CEN03000 

C199 SUM=SUM+ ( 1 . / ( (X**4) -0 . 25) ) CEN03010 

C EXACT=SUM*D 1 CEN03020 

CEN03030 

L I N E =3 CEN03040 

I F (I COUNT. EQ. 50) GO TO 1 4 1 CEN03050 

GO TO 143 CEN03060 

141 WRITE (2, 175)TIME,U(N/2) .LINE CEN03070 

C D I F= (U (N/2) -EXACT) /EXACT CEN03080 

C IF (DABS(DIF) .LE.0.15)WRITE (2,177) CEN03090 

I C0UNT=O CEN03100 

CEN031 10 

143 DO 150 1 = 1 , N CEN03120 

UTN ( I ) =UT ( I ) CEN03130 

150 UT(I)-U(I) CEN03140 

CEN03150 

IF (TIME .GT.TT) GO TO 500 CEN03160 

CEN03170 

GO TO 130 CEN03180 

175 FORMAT (F10.8, 1X.F10.8, IX, 11) CEN03190 

176 FORMAT (/IX, 1 TIME DEFLECTI ON AT L/2 M5X, ' EXACT' ) CEN03200 

177 FORMAT (' — CEN03210 

6 ') CEN03220 

500 STOP CEN03230 

END ' CEN03240 

CEN03250 

CEN03260 

SUBROUTINE I NVERT (AO , A, N) CEN03270 

DOUBLE PRECISION A (70, 70) , AO (70, 70) CEN03280 

CEN03290 

DO 1 1=1, N CEN03300 
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DO 1 J-l.N 

CEN03310 

1 

A (1 , J) =A0 (1 , J) 

CEN03320 

CEN03330 


NP=N+1 

CEN03340 


A ( 1 , NP) =1.0 

CEN03350 


DO 10 1=2, N 

CEN03360 

10 

A (1 , NP) =0.0 

CEN03370 

CEN03380 


DO 40 J= 1 , N 

CEN03390 


DO 20 LX=1 ,N 

CEN03400 

20 

A (NP , LX) =A (1 , LX+1) / A ( 1 , 1) 

CEN03410 


DO 30 KX=2 , N 

CEN03420 


DO 30 LX=1 , N 

CEN03430 

30 

A (KX — 1 , LX) -A (KX, LX+1) -A (KX, 1 ) *A (NP,LX) 

CEN03440 


DO 40 L X= 1 , N 

CEN03450 

40 

A (N , LX) =A (NP , LX) 

CEN03460 

CEN03470 


RETURN 

CEN03480 


END 

CEN03490 

CEN03500 


SUBROUTINE JACOBI (K,M,N, 1 FPR, X, E 1 GV) 

CEN03510 

C 

SUBROUTINE JACOBI 

CEN03520 


IMPLICIT REAL*8(A-H,0-Z) 

CEN03530 


DOUBLE PRECISION A (70, 70) , B (70, 70) , X (70 , 70) , E 1 GV (70) , D (70) , 

CEN03540 

6K(70,70) ,M (70, 70) 

CEN03550 


IFPR=0 

CENO 356 O 

C 

COMMON/K.M/ 

CEN03570 

C 

WRITE (2, 1051) 

CENO 358 O 

ci 051 

FORMAT (/IX , 1 INPUT DATA ') 

CEN03590 

C 

READ (1 ,*) N, 1 FPR 

CENO 36 OO 

C 

WRITE (2, 1001) N, 1 FPR 

CENO 36 IO 

C 

DO 1010 1 = 1 , N 

CEN03620 

C 

READ (1 ,*) (A (1 , J) ,J=1 ,N) 

CEN03630 

C 

WRITE (2, 1 1 10) (A(l ,J) ,J=1 ,N) 

CEN03640 

C1010 

CONTINUE 

CENO 365 O 

c 

DO 1020 1 = 1 , N 

CENO 366 O 

c 

READ (1 ,*) (B (1 ,J) ,J=1 ,N) 

CENO 367 O 

c 

WRITE (2, 1 110) (B ( 1 , J) ,J=1 , N) 

CENO 368 O 

Cl 020 

CONTINUE 

CEN03690 

Cl 001 

FORMAT (21 10) 

CEN03700 

ci no 

FORMAT (8F 10.4) 

CEN03710 


DO 2 1=1, N 

CEN03720 


DO 1 J-l.N 

CEN03730 


A ( 1 , J) =K 

CEN03740 


B ( 1 , J) =M ( 1 , J) 

CEN03750 

1 

CONTINUE 

CEN03760 

2 

CONTINUE 

CEN03770 


NSMAX=1 5 

CEN03780 

c 

WRITE (2, 1980) 

CEN03790 

1980 

FORMAT (/IX, 1 EIGENVALUES ') 

CENO 38 OO 


RT0L=1 . D- 1 2 

CEN03810 


1 0UT=2 

CEN03820 


DO 10 1=1, N 

CEN03830 

c 

PRINT*, 'FLAG ',1,' A = ' , A ( 1 . 1 ) , ' B = 1 . B ( 1 , 1 ) 

CEN03840 


IF (A (1 , 1) .GT.O.AND.B (1 , 1) .GT.O.) GO TO 4 , 

CEN03850 
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WRITE (I OUT, 2020) 

STOP 

4 D(I)=A(I,I)/B(I,I) 

10 E IGV (I) =D (I) 

DO 30 1=1, N 
DO 20 J=1 ,N 
20 X ( I , J) =0 . 

30 X(l,l) = 1.0 

IF (N.EQ.1) RETURN 
C 

C INITIALIZE SWEEP COUNTER AND EIGEN ITERATION 
C 

NSWEEP=0 
NR=N- 1 

40 NSWEEP=NSWEEP+1 

I F (IFPR.EQ. 1) WRITE ( I OUT, 2000) NSWEEP 
PRINT*,' SWEEP NUMBER... '.NSWEEP 
C 

C CHECK IF PRESENT OFF DIAGONAL ELEMENT IS TOO LARGE 
C 

EPS= (0.01 **NSWEEP) **2 
DO 210 J=1,NR 
JJ=J+1 

DO 210 K1=JJ,N 

IF (DABS (A (J,K1)) . LT . 1 . D-20) GO TO 211 
EPTOLA® (A (J , K 1) *A (J , K 1 ) ) / (A (J , J) *A (K 1 , K 1) ) 

GO TO 212 

211 EPTOLA=0.0 

2 1 2 EPTOLB= (B (J , K 1 ) *B (J , K 1 ) ) / (B (J , J) *B (K 1 , K 1 ) ) 

IF ((EPTOLA.LT.EPS) .AND. (EPTOLB.LT. EPS) ) GO TO 210 
AKK=A (K 1 , K 1 ) *B (J , K 1 ) -B (K 1 , K 1 ) *A ( J , K 1 ) 

A J J=A (J , J) *B (J , K 1 ) -B (J , J) *A (J , K 1 ) 

AB=A (J, J) *B (K 1 ,K1) -A (K 1 , K 1 ) *B (J , J) 

CHECK= (AB*AB+4.*AKK*AJJ) /4. 

C PRINT*, 'THIS IS CHECK... '.CHECK 
I F (CHECK) 50,60,60 
50 WRITE (I OUT, 2020) 

STOP 

60 SQCH=DSQRT (CHECK) 

D l=AB/2 .+SQCH 
D2=AB/2 . -SQCH 
D EN=D 1 

I F (DABS (D2) .GT.DABS (Dl) ) DEN=D2 
I F (DEN) 80,70,80 
70 CA=0 . 

CG= (-1 .) *A (J , K 1) /A (K 1 , K 1) 

GO TO 90 
80 CA=AKK/DEN 

CG= (- 1 .) *AJ J/DEN 
90 I F (N-2) 100, 190, 100 
100 JP1=J+1 

JM1 =J- 1 
KP1=K1+1 
KM1 =K 1 - 1 

IF(JMl-l) 130,110,110 


CEN03860 
CEN03870 
CEN03880 
CEN03890 
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110 DO 120 1=1, JM1 
AJ=A (I , J) 

BJ=B ( I ,J) 

AK=A (1 , Kl) 

BK=B (I ,K1) 

A (I ,J)=AJ+CG*AK 
B (I , J) =BJ+CG*BK 
A ( I , K 1) =AK+CA*AJ 
120 B ( I , K 1) =BK+CA*BJ 
130 I F (KP1-N) 140, 140, 160 
140 DO 150 I =KP 1 , N 
AJ=A (J , I) 

BJ=B (J , I ) 

AK=A (K 1,1) 

BK=B (K 1 , I ) 

A (J, I) =AJ+CG*AK 
B (J , I ) =BJ+CG*BK 
A (Kl , I) =AK+CA*AJ 
150 B (K 1 , I ) =BK+CA*BJ 
160 IF (JP1-KM1) 170,170,190 
170 DO 180 I=JP1 ,KM1 
AJ=A (J, I) 

BJ=B (J, I) 

AK=A (I , K 1) 

BK=B (I ,K1) 

A (J, I ) =AJ+CG*AK 
B (J , I ) =BJ+CG*BK 
A (I , Kl) =AK+CA*AJ 
180 B (I , Kl) =BK+CA*BJ 
190 AK=A (Kl ,K1) 

BK=B (Kl ,K1) 

A (K 1 , K 1 ) =AK+2 . *CA*A ( J , K 1 ) +CA*CA*A ( J , J) 

B (Kl , K 1 ) =BK+2 . *CA*B (J , Kl) +CA*CA*B (J , J) 

A (J , J) =A (J , J) +2 . *CG*A (J,K1)+CG*CG*AK 
B (J , J) =B (J , J) +2 . *CG*B (J , K 1) +CG*CG*BK 
A (J , Kl) =0 . 

B (J , K 1) =0 . 

UPDATE EIGENVECTOR MATRIX 

DO 200 1=1, N 
XJ=X (I , J) 

XK=X (I , K 1) 

X(l ,J)=XJ+CG*XK 
200 X (I , Kl) =XK+CA*XJ 
210 CONTINUE 
C 

C UPDATE EIGENVALUES 
C 

DO 220 1=1, N 

I F (A (I , I) .GT.O.AND.B (I , I) .GT.O) GO TO 220 
WRITE (I OUT, 2020) 

STOP 

220 E I GV ( I ) =A (1,1) (1,1) 

I F (I FPR.EQ.O) GO TO 230 
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WRITE (I OUT, 2030) 

WRITE (I OUT, 20 10) (E I GV (I ) , I = 1 ,N) 

C 

C CHECK FOR CONVERGENCE 
C 

230 DO 240 1-1, N 
TOL=RTOL*D ( I ) 

D I F=DABS (EIGV(I)-D(I)) 

IF (DIF.GT.TOL)GO TO 280 
240 CONTINUE 
C 

C CHECK ALL OFF D I AG ELEMENTS TO SEE IF ANOTHER SWEEP IS REQ'D 
C PRINT*,' RTOL 1 ,RTOL 

EPS=RT0L**2 
DO 250 J= 1 , NR 
JJ=J+1 

DO 250 K1=JJ,N 

I F (DABS (A (J,K1)) .LT. 1 .D-30)GO TO 251 
EPSA= (A (J , K 1 ) *A (J , K 1 ) ) / (A (J , J) *A (K 1 , K 1 ) ) 

GO TO 252 

251 EPSA=0.0 

C PRINT*, 1 EPSA '.EPSA,' EPS '.EPS 

252 EPSB= (B (J , K 1 ) *B (J , K 1 ) ) / (B (J , J) *B (K 1 , K 1 ) ) 

C PRINT*,' EPSB ' , EPSB, 1 EPS '.EPS 

IF ( (EPSA. LT. EPS) - AND . (EPSB.LT.EPS) ) GO TO 250 
GO TO 280 
250 CONTINUE 
C 

C FILL OUT BOTTOM TRIANGLE OF RESULTANT MATRICES 6 SCALE EIGENVECTORS 
C 

255 DO 260 1=1, N 
DO 260 J=1 ,N 
A (J , I ) =A ( I , J) 

260 B (J , I ) =B ( I , J) 

DO 270 J=1 , N 
BB=DSQRT (B (J , J) ) 

DO 270 Kl-l.N 
270 X (K 1 , J) =X (K 1 , J) /BB 
C WRITE (IOUT, 310) 

C DO 300 1=1, N 

C300 WRITE (IOUT, 2010) (X (I ,J) ,J=1 ,N) 

310 FORMAT (/IX, ' THE EIGENVECTORS ARE ') 


RETURN 

UPDATE THE ’ D ' MATRIX AND START NEW SWEEP IF ALLOWED 

280 DO 290 1=1, N 
290 D ( I ) =E I GV ( 1 ) 

IF (NSWEEP.LT. NSMAX) GO TO 40 
GO TO 255 

2000 FORMAT (/IX,' SWEEP NUMBER IN JACOBI = ' , 1 4) 
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2010 FORMAT (/IX, 6E20. 12) 

2020 FORMAT (/IX,' **** ERROR SOLUTION STOP / MATRICES NOT POSITIVE 

6DEFINITE') 

2030 FORMAT (/IX, 1 CURRENT EIGENVALUES IN JACOBI ARE ') 

END 


SUBROUTINE DAMP (N , E I GV, X ,M, DAMRAT , C) 

IMPLICIT REAL*8(A-H,0-Z) 

DOUBLE PRECISION X (70 , 70) ,T (70, 70) , M (70 , 70) , C (70, 70) , E 
SAT (70) 

DO 10 1=1, N 
E I GV ( I ) =DSQRT (E I GV ( I ) ) 

DO 10 J=1,N 
10 C ( I , J) =0 .0 

DO 20 I 1=1 , N 

DA=2 . *DAMRAT (I I ) *E I GV ( I I ) 

DO 20 1-1 .N 
DO 20 J=1,N 

20 C (I ,J)=C(I ,J)+X(I , I l)*X(J, I l)*DA 

DO 30 1-1, N 
DO 30 J=1,N 
T (I , J) =0.0 
DO 30 Kl-l.N 

T ( I , J) =T ( I , J) +M ( I , K 1 ) *C (K 1 , J) 

DO 40 1=1, N 
DO 40 J=1 ,N 
C ( I , J) =0 .0 
DO 40 K 1 = 1 ,N 

40 C(l , J) =C (I ,J)+T(I ,K1)*M(K1,J) 

C DO 50 1 = 1 » N 

C50 WRITE (2, 120) (C (I ,J) , J=1 ,N) 

120 FORMAT (6D 14. 4) 

RETURN 

END 


30 
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THIS IS CENTRAL DIFFERENCE SOLUTION 


STIFFNESS — 0.531100006D+02 

DAMPING 0. 130999982D-01 

FORCE 0.400000066D-02 

FREQUENCY — 0 . 251 300049D+02 

TIME DEFLECTION AT L/2 


0.00500000 0.07560742 1 
0.01000000 0.07411487 1 
0.01500000 0.07822479 1 
0.02000000 0.07755833 1 
0.02500000 0.07321689 1 
0.03000000 0.07439081 1 
0.03500001 0.07596175 1 
0.04000001 0.07030083 1 
0.04500001 0.06480695 1 
0.05000001 0.06403831 1 
0.05500001 O.05838368 1 
0.06000001 0.04628549 1 
O.O65OOOOI 0.03770042 1 
0.07000001 0.03044851 1 
0.07500001 0.01549482 1 
0.08000001 -0.00155219 1 
O.O85OOOOI -0.01381828 1 
0.09000001 -0.02866951 1 
0.09500002 -0.04917751 1 
0.10000002 -0.06667677 1 
0.10500002 -0.08058453 1 
0.11000002 -O.0980749I 1 
0.11500002 -O.II6OIO54 1 

0.12000002 -0.12788046 1 
0.12500002 -0.13797468 l 
0.13000002 -0.14964225 1 
0.13500002 -0.15657816 1 
0.14000002 -0.15730855 1 
0.14500002 -0.15770900 1 
0.15000002 -0.15621028 1 
0.15500003 -0.14753221 1 
0.16000003 -0.13531903 1 
0.16500003 -0.12327482 1 
0.17000003 -0.10672159 1 
0.17500003 -0.08461688 1 
0.18000003 -0.06275240 l 
0.18500003 -0.04085973 1 
0.19000003 -0.01462197 1 
0.19500003 0.01269540 1 
0.20000003 0.03673648 1 
0.20500003 0.06085743 1 
0.21000003 0.08652612 1 
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